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In this paper we study nonequilibrium dynamics of one dimensional Bose gas from the general 
perspective of dynamics of integrable systems. After outlining and critically reviewing methods 
based on inverse scattering transform, intertwining operators, q-deformed objects, and extended 
dynamical conformal symmetry, we focus on the form-factor based approach. Motivated by possible 
applications in nonlinear quantum optics and experiments with ultracold atoms, we concentrate 
on the regime of strong repulsive interactions. We consider dynamical evolution starting from 
two initial states: a condensate of particles in a state with zero momentum and a condensate of 
particles in a gaussian wavepacket in real space. Combining the form-factor approach with the 
method of intertwining operator we develop a numerical procedure which allows explicit summation 
over intermediate states and analysis of the time evolution of non-local density-density correlation 
functions. In both cases we observe a tendency toward formation of crystal-like correlations at 
intermediate time scales. 

I. INTRODUCTION 

Quantum many-body systems in low dimensions can not typically be described using mean-field approaches. This 
makes analysis of their non-equilibrium dynamics particularly challenging. However following demand from recent 
experiments, certain progress has been achieved in developing theoretical methods which can address this problem. In 
this paper we use Bethe ansatz solution to study nonequilibrium dynamics of the one dimensional Bose gas interacting 
via contact interaction. This microscopic model of i5-interacting bosons is called the Lieb-Liniger (LL) model 
It belongs to a general class of models which can be studied using Bethe ansatz solution [1] . These exactly solvable 
models are characterized by an infinite number (in the thermodynamic limit) of conserved quantities which originate 
from the nature of collision processes. In the LL model the two particle collisions can not change momenta of scattering 
particles but only give rise to a phase shift. Moreover any many-body collision is factorizable into a sequence of two- 
body scattering events. While implications of the exact solution for thermodynamic properties have been discussed 
before [U, HQ, their consequences for non-equilibrium dynamics are not known. This will be the central question of 
our paper. 

Systems of one dimensional Bose gases with contact interactions have been recently realized using ultracold atoms 
[Toj and current experiments allow wide control over parameters of the microscopic Hamiltonian (see [ll[ , for recent 
review). For example, the effective strength of the repulsive interaction can be tuned either by changing the density 
of the atomic cloud, or by modifying the strength of the transverse confinement, or by changing the scattering length 
using magnetically tuned Feshbach resonance. In particular a very interesting regime can be achieved when repulsion is 
so strong that bosonic atoms become essentially impenetrable 13 1 and the system undergoes effective fermionizationf^ 



|8| . Recent studies of thermodynamic properties demonstrated excellent agreement between experiments and the exact 
solution [1, . However parameters of the Hamiltonian can also be changed dynamically. In this case one needs to 
study dynamics starting with the initial state which is not an eigenstate of the Hamiltonian. Generally such dynamics 
• y—{ involves superposition of coherent evolution of all eigenstates of the system. This is the problem that we address in 
this paper. 

^ Another class of systems where the LL model appears naturally is light propagation in one dimensional photonic 

- - ■ fibers. It has been known since the sixties (see e.g. ref [l3|) that propagation of classical pulses of electromagnetic 
waves in one dimensional nonlinear medium can be described by the nonlinear Schroedinger equation, which can be 
interpreted as the operator equation of motion arising from the LL microscopic Hamiltonian. In this case nonlinearity 
is proportional to the nonlinear polarizability of the medium, x^*^^ (see e.g. Ref [l5?]). Typical propagation problem in 
one dimensional nonlinear fiber corresponds to the classical limit of the quantum LL model subject to nonequilibrium 
boundary and/or initial conditions [16]. Earlier analysis assumed that photon nonlincarities come from interaction of 
photons with two level atoms, in which case nonlinearity corresponds to the attractive interaction in the LL model [l7| . 
Corresponding quantum models have bound states which eventually lead to formation of solitons and classical regime 
[l^ . Earlier analysis of optical systems also assumed regime of weak nonlincarities since increasing nonlincarities 
in two level systems would lead to strong photon losses. However recent techniques utilizing electromagnetically 
induced transparency [T^ make it possible not only to realize effective repulsion between photons but also to increase 
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dramatically the strength of nonlinearities(20j. Optical realizations of the LL model should allow direct measurements 
of both the first and the second order coherences [21[, g^^\t, r) = {"^'^t)"^!{t + T)) and g^'^\t, r) — {p(t) p{t + t)) , where 
'i'{t) is the amplitude of the propagating electromagnetic field and p{t) — is the field intensity, measuring 

the number of photons. 

There are three main reasons why we undertake a detailed analysis of the non-equilibrium quantum dynamics of 
the LL model. The first reason is that dynamics of integrable models, such as the LL model, should exhibit special 
features originating from an infinite number of integrals of motion. One manifestation of constraints on the phase 
space available for dynamics should be the absence of thermalization and ergodicity ^ ,22,. .23,] . Conservation laws 
originating from integrability also lead to the so-called Mazur inequalities [l^l , which should result in special transport 
properties. For example, the possibility of an infinite Drude weight has been discussed for several integrable models 
of lattice fermions [25|. Our second motivation for analysis of the LL model is to use it as a testing ground for 
developing new theoretical methods and techniques, which can then be applied to a broader class of models and 
systems. Our work should provide a general framework for understanding dynamics of integrable systems and for 
future developments of non-perturbative methods in the study of non-equilibrium dynamics. We emphasize that our 
analysis uses special properties of exactly solvable models, and the methods we apply here are fundamentally different 
from the more conventional perturbative approaches discussed in Refs. (26j.[27j. Finally we point out that recent 
progress in the areas of ultracold atoms, quantum optics, and low-dimensional strongly-correlated materials makes it 
possible to fabricate concrete physical systems, which can be accurately described by the LL model. Thus our third 
reason for analyzing this model is that theoretical predictions for the time-dependent evolution of correlation functions 
can be measured directly in experiments. Specific system in which one can realize dynamical experiments discussed in 
our paper are one dimensional Bose gases in optical lattices and magnetic microtraps (see Refs. [l^H^ for a review). 
Realization of such experiments should provide a unique opportunity to study non-equilibrium dynamics of strongly 
correlated exactly solvable systems. 

This paper is organized as follows. Section II provides a brief critical overview of several different methods which 
can be used to describe non-equilibrium dynamics of integrable quantum systems such as the LL model. Section 
III provides an in-depth discussion of one of these approaches, the so-called form-factor technique, which relies on 
numerical summation over intermediate states. The overlap with the initial states in the strongly interacting regime 
is found using the method of intertwining operator. In section IV we focus on the non-equilibrium dynamics of the 
many-body system described by the LL model. We calculate the time evolution of the density-density correlation 
function for two different condensate-like initial states: when all particles are in a state with zero momentum and 
when all particles are in a Gaussian wavepacket in real space. Summary of our results and conclusions are given in 
Section V. To provide additional background to the readers, in appendices we provide basic facts about the algebraic 
Bethe Ansatz formalism and the inverse scattering transform for the model of one dimensional Bose gas with both 
repulsive and attractive interactions. 

II. A REVIEW OF METHODS TO INVESTIGATE NONEQUILIBRIUM DYNAMICS OF 

INTEGRABLE MODELS 

In this paper we analyze the time evolution of correlation functions in integrable models when the initial state is 
not an eigenstate of the Hamiltonian. Here we give an overview of (some of) the possible approaches to studying 
nonequilibrium dynamics of integrable systems. We only discuss approaches which explicitely use the property of 
integrability and exclude more conventional techniques which rely on approximations and perturbative expansions, 
such as the Keldysh formalism . Although the focus of this paper is on the ID Bose gas, methods discussed here 
can be applied to other exactly-solvable models, e.g. for the Gaudin-type integrable models 

The Hamiltonian of the ID interacting Bose gas in a finite size system has the following form Q 



Here ^(x) is a bosonic field in the second quantized notations, c is an interaction constant, L is the size of the system. 
The equation of motion for this Hamiltonian is the Nonlinear Schrodinger (NS) equation. The first quantized version 
of this Hamiltonian corresponds to the Lieb-Liniger model of Bose gas interacting with the contact interaction. The 
Bethe ansatz solution of this model has been given in refs. [H,[l|i[lli[l|- 

Below we discuss several possible approaches to the study of nonequilibrium dynamics of nonlinear integrable 
models, such as the LL model ([T]). 




(1) 
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A. Quantum inverse scattering method 

One approach to analyzing dynamics of model ([I} would be to use the formalism of the inverse scattering problem. 
This method relies on the solution of the quantum version of the Gelfand-Levitan-Marchenko equation (see Appendix 



A). From the Inverse Scattering Transform (see e.g. |3l| for review) it is known that the solution is given by the 
following infinite series representation 



N , N 



*(^) = E / n ^ n -^SNm. {k}; x)R\p^) ■ ■ ■ R\p^)R{k^) ■ ■ ■ R{ko) (2) 

N=0'' 1=1 j=0 ^ 

where the operators R{p) diagonalize the problem on the infinite interval [—oo, oo]. The function can be cast into 
different forms [32] . One specific representation is given by 

9n{{p}, {k}; x) = — ^- — (3) 

The above perturbative expansion is a quantum version of Rosales expansion [33j.[3^. Thus, 

The inverse expression (direct Gelfand-Levitan transform) is also available (see [35l| for the finite interval), 

i?t(^)===^= [ dx-i'^x)e-"'^ + ^= [ da;ida;2rfa;352(a;i,a;2,X3;g)*^(a;i)*t(a;2)^f(a;3) + ... (5) 
VZTT J yzn J 

where e.g. ^ g2{xi,X2,X3-q) = c9{x2 - X3)9{x3 - xi) exp{-iq{xi + X2 - X3)), g4{xi,X2,X3,X4^,X5;q) = c^Oixs - 
X5)0{x5 — X2)9{x2 — X4)0{x4 — xi) e-K.p{-'iq{xi + X2 + X3 — X4 — 2:5)). Here the operators of refiection coefficient 
i?(^),i?^(f) diagonalize the Hamiltonian and satisfy the Zamolodchikov-Faddeev algebra (see Appendix A for more 
details) 

[H,R^{0] = eRHO, (6) 

Rmio = s{e - ORiom), Rm\a = ^(c - mHORio + '^^^{i - o- (7) 

where the scattering matrix is 5'(^ — ^') = {£, — £,' — ic) — S,' + ic) . These relations make the problem of finding time 
evolution easy: the factor exp[zx(EQ fci — X^i Pi)] i'^ Eq. ([3]) simply needs to be replaced by the factor exp[ix(X]o ^i — 
Pi) ~ *^(So^ ~ -Pf)]- Application of this formalism to non-equilibrium problems is also straightforward: 
one has to decompose the initial state, written in terms of bosonic \I'(a;)-operators, into a series of Zamolodchikov 
-Faddeev operators R{^), generated by the inverse transform, and then find the time evolution according to ([6]) using 
the direct transform. 

Strictly speaking, this expansion is valid only for an infinite interval. For finite intervals it gives rise to singularities 
in expressions for the correlation functions, which, however, can be corrected by careful consideration of expressions 
in each order [3^. This series can be summed up explicitely only for the infinite value of the coupling c [37|.[3^. 
Unfortunately this approach is very difficult to use for calculating time dependence of the correlation functions at 
finite c even in equilibrium. We are not aware of the reproducibility of the asymptotic results which would correspond 
to the Luttinger liquid power-laws. On the other hand the advantage of this approach is the possibility to generalize 
it to other specific boundary or initial conditions. In principle this formalism should allow to include impurities [ssj ] 
and can be extended to multi-component generalizations of the NSE [39j . 



B. Intertwining operator 

Some integrable (and sometime quasi-exactly solvable) models have the following property 

IHo = HJ (8) 

where two different Hamiltonians 11^ ^ are connected (intertwined) by the action of some operator /, called the 
intertwining operator. In general, if -ffo and He belong to two different representations of the same algebra, then 
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intertwining operators exchange these two representations of the same algebras (or, saying mathematically, establish 
a certain homomorphism between them). If Hq is a Hamiltonian for free, noninteracting particles, and He is the 
Hamiltonian for the interacting system with interaction constant given by c, then the evolution of observables in 
the interacting model can be related to the evolution of the non-interacting one. Hence dynamics of the interacting 
Hamiltonian can be mapped to dynamics of the noninteracting Hamiltonian using 

e'"-' ^ le'^^^rK (9) 

Few comments are in order. The existence of the intertwining operator for integrable models follows from the 
following facts: (i) the Bethe states of all integrable models are parametrized by the integer numbers which label the 
eigenfunctions of noninteracting problem. Therefore the wavenumbers of interacting model are analytically connected 
to the wavenumbers of the noninteracting one; (ii) in the language of the coordinate Bethe ansatz a system of 
N particles is described by the wave function defined in the A^-dimensional space divided by the hyperplanes on 
which the collision processes occur. Outside these hyperplanes a system behaves as noninteracting. Transition 
amplitudes between A^! different regions (outer space of the hyperplanes) are the same for all eigenstates because of 
the permutation symmetry. These transition amplitudes define a unitary transformation which is nothing but the 
intertwining operator. However, one should note that in general the operator / is not unitary. Thus, the wavefunctions 
of exactly solvable models are not orthonormalized, their overlaps are given by the determinants of some matrix via 
Gaudin-Slavnov formula [|,[3. Defining orthonormalized wavefunctions one can construct a unitary version of the 
operator /. 



1. Degenerate ajfine Hecke algebra: first quantized notations 

In the case of the Lieb-Liniger model the intertwining operator can be related to the representation of the degenerate 
affine Hecke algebra (ioj . For the XXZ spin chains there should exist similar intertwining operator between different 
representations of the Uqi^sh) and Temperley-Lieb algebras. One reason to expect this is because both LL model and 
XXX spin chain share essentially the same i?-matrix satisfying the Yang-Baxter equation. Some earlier discussion on 
the XXZ chain are contained in 41]. 

The construction of / for the ID Bose gas relies on the notion of Dunkl operator 

d, = -id., + i^ ^(e(x, - Xj) - l)sij + i^ ^(e(xi - xj) + l)s,,j (10) 

j<i j>i 

where e(a:) is a signature function, operator Si,j provides a representation of the Artin's relations for the braids 
and exchanges coordinates of particles i and j and satisfy The operators di,Sj = Sjj+i form a 

representation of the degenerate affine Hecke algebra. Another representation of the same algebra is formed by the 
ordinary differential (difference) operator —idi and the integral operator Qi [4^ representing scattering matrix and 
acting on the arbitrary function /(..., Xi,Xi+i, . . .) as 

Qif{. . .,Xi,Xi+i, ...)=/(.. . ...)- c f{...,Xi- t,Xi+i +t,.. .)dt (11) 

Jo 

The intertwining operator / interchanges these two representations of the affine Hecke algebra, (c?i, Si^j) and {—idi, Qi) 
and thus intertwines the Dunkl operator di and the ordinary partial differential operator. Explicitly 

^ e{x^-l(l) < . . .X^-l(N))s^-lQw (12) 

where s^-i = si^ . . . Sj^Sj , an d Qw = QiiQi2 ■ ■ ■ Qip (1 < ii, • ■ • , ip < ~ 1) and where w is a transposition from 
the symmetric group Sn |40| . 

All conserved quantities X„ for the ID Bose gas system are given by the powers of the Dunkl operator, I„ — J^i '^{dl'') 
where 7r(.) is projection onto symmetric subspace. Thus, the Hamiltonian ^ is just equal to X2. 

As a side remark we note that these facts are convenient for the formulation of the generalized Gibbs ensemble 
(GGE) approach to non-equilibrium dynamics of the LL model [l^]. In particular, the GGE density matrix of the 
ID Bose gas should be related to the noninteracting density matrix by the intertwining relation. As discussed in 
Ref. [4^, the GGE conjecture should work when eigenvalues of the integrals of motion are parametrized by either 
{0, 1} (fermionic-like systems) or by integers {0, 1,2,.. .} (bosonic-like systems). The case of the Bose gas belongs to 
the second class. We therefore conjecture that the GGE is applicable to the ID Bose gas (LL model), at least for the 
local observables. 
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2. Graphical representation: Gutkin approach 

An interesting approach to the ID Bose gas based on the intertwining operator has been developed by E. Gutkin in 
a series of papers [4l| , [12] , [11] ; @ ■ 1^ ^^'^ later applied to the problem of propagation of an optical pulse in nonlinear 
media in Ref. [46| . 

In this approach the intertwining operator is expanded as a series of different contributions labeled by the so-called 
collision graphs [93|- The evolution of the quantum field '^{x,t) governed by the NSE is given entirely in terms of 
the evolution of the field ^o{x,t) for the free Schrodinger equation, idt'^Q{x,t) = —dl'^o{x,t), and by the quantities 
Or (defined in terms of some distributions) entering the expansion of the intertwining operator over a set of collision 
graphs 

1= d'^x d''yarixi,X2r-- ,Xn,yl,y2r■■yn)'^Uxl)^lix2)■■■'^'lixr^)'^o{yl)'^ (13) 

r,n=g(r)"^ 

and similarly for the operator with ar replaced by the related distribution br- There is a special choice of ar and 
br for which the operator / is unitary. For details about this approach we refer to the original papers by Gutkin. 

The advantage of this approach is the possibility to find the propagator of the nonlinear problem and thus, poten- 
tially, solve the initial state evolution problem for an arbitrary initial state. Computationally one only needs to deal 
with the free fields ^'o(a^) and their time evolution with a free (noninteracting) Hamiltonian. We note that the collision 
graph expansion is highly non-perturbative and can not be rewritten as an expansion in powers of c in opposition to 
the Gelfand-Levitan-based expansion of Sec. Ill Al A disadvantage of this approach is rapidly growing complexity of 
collision graph expansion with increasing number of particles. Also it is not known at present whether some of the 
most relevant graphs can be summed up in order to develop useful approximation schemes. 



3. Second quantized approach 

In series of papers [131 Sasaki and Kebukawa developed a field-theoretical approach to the LL model. Although 
they did not present it in this context, the second quantized form of the intertwining operator can be recognized in 
their construction fOl]. In this approach one considers the LL Hamiltonian written in the usual second-quantized 
form, 

2 ^ 

^^ = Y1 2^°p"P + ^'^l+ral-raqap. (14) 

Here aL Oj, are bosonic creation and annihilation operators corresponding to momenta pi — 2nfmi/L {rii should be 
integer). The second-quantized analog of the Bethe ansatz wavefunction for N particles has the form 

N 

|*gi,g2,...9iv> = ^gi,g2,...9iv 11 ^^^^ ^ ' ^^'^'^ IT °E"=i ■ . P. , +9.'^'^ ^^^^ 

Pij;l<i<j<N l<i<j<N i=l ' 

where Bq-^^q^ ^ ^q^^ is a normalization factor, pij ~ Pi—Pj = '^Pj.i and (i = 1, . . . , N) are defined as an integerx27r?i/L, 
d(Pi,j] ki_j) — —ki_j/(pi,j — ki_j) and kij — —kj^i are solutions of the Bethe ansatz equations 

cot{^^) ^—{h- kj) = —{2h,, + V (fc,,, - fcj- i) + - g,), l<i<3<N. (16) 

The eigenvalues in these notations are given by Eq^^,,,^q^ — X^iLi k'f/2m = Y^^=i J2j^ii^i,j + '?i)^/2m. 

It is interesting that authors of [J?] found a unitary operator which transforms eigenstates of noninteracting system 
Hq into eigenstates of the interacting Hamiltonian He- In the second quantized form this unitary operator is given by 

n n 

Y Bq,^q,^,„^q^Aq,^q,^,„^q^ [| ^(^'^J ^ '^^.J ) ll "E?=i +9, ll ^^'^^ 

gi <(/2 ■ ■ -^Ti Pi. j ; 1^^ *<Jf^^^ l<i<j<n i—1 i—1 

where are normalization factors for non-interacting eigenstates. Explicit calculations show that this unitary 

operator is indeed equivalent to the intertwining operator from previous subsections. 

This approach has advantages for initial states which can be readily represented using formalism of second quanti- 
zation. For the case of a large coupling constant it gives the same result as the one used later on in the text. 



6 



C. q-bosons 

In a series of recent papers [11], [i^ a system of g-bosons hopping on ID lattice has been studied. This model is 
defined by its Hamiltonian 



1 

-^(&t6„+i+5„&l+i-27V„) (18) 



where the periodic boundary condition is imposed. Operators i3„, B^^ and A^„ satisfy the g-boson algebra 

[N„b]] = blS,,, [N,,bj]=-bAj, [b,,b]] =q-^^^5,j (19) 

where q — e"' . The operator of the total number of particles TV — X^jli commutes with the Hamiltonian Hq. 
The continuum limit of the model is defined using the limiting procedure 

cS 

When applied to the system of (/-bosons this procedure gives 

H=[ dx [d^b'' {x)d^b{x) + cb^ {x)b\x)b{x)bix)] , (21) 
Jo 

where [b{x),b^y)] = 6{x — y) are the canonical Bosc fields. Therefore ID Bose gas with contact interaction can be 
regarded as a continuum limit of the g-boson lattice model. Apparently a solution of the inverse scattering problem 
for the NSE should be related to the one for the (/-boson model. An explicit form of this relationship is not known, 
although should exist '92] . 

This approach can be used to study non-equilibrium dynamics. Solutions of particular evolution problem for q- 
bosons can be transformed into solutions of NS problem via the limiting procedure (|20|) . An attempt to construct 
the evolution operator for a single-mode g-boson system has been made in [s^]. Even for this simple case the g-path 
integral has an involved structure which includes integration with the measure corresponding to non-flat phase space 
and complicated action having non-trivial dynamical phase. All this makes the possibility of explicit evaluation of 
the g-path integral questionable. It is interesting to note that dynamics in this phase space is naturally constrained 
by integrals of motion. This is reminiscent of recent suggestions of the importance of GGE for dynamics of integrable 
models [H. 



D. Extended conformal symmetry 



It is known that the low-energy description of the ID Bose gas can be formulated as a bosonic Luttinger liquid 
[Sll [52! . This description operates with a linear spectrum of bosonic modes and treats the interaction part essentially 
exactly. From the more formal point of view, in this low-energy description the thermodynamic limit of the NS 
system represents a conformal field theory (GET) with the central charge c = 1. This means, in particular, that the 
low-energy (bosonic Luttinger) Hamiltonian is a linear combination of the zero-momentum generators Lqj ^0 (for the 
right and for the left moving parts) of the conformal symmetry algebra, i.e. the Virasoro algebra 

[Ln,L„,] = [n - m)Ln+m + Y^"("^ ^ l)(5„+m,o (22) 

and the same for Z„, ([i„, Lm] — 0, Vn, m). The space of states is the so-called Verma module (representation space) 
for this Virasoro algebra. The Kac table defines a conformal tower for this c = 1 GET which allows to compute any 
correlation functions. One could also include finite size corrections and effects of irrelevant perturbations ^3], 54]. 

The non-equilibrium dynamics can not be correctly described in terms of the low-energy, Luttinger liquid theory 
only. The reason is two-fold: first, if the initial state represents a "large" perturbation over equilibrium, which includes 
the overlap with the whole spectrum, the excitations involving nonlinear part of the spectrum should be important. 
Second, even for relatively weak perturbations, which do not thermalize because of integrability, the dynamics at large 
times leads to growing correlations over the entire momentum space. So irrelevant contributions from the interaction 
part of the Hamiltonian should become important at long times. 
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It appears that nonlinearity of the dispersion as weh as irrelevant (from the equihbrium point of view) parts of 
interaction can be included using the extended conformal symmetry called Wi+oo- The Wi+oo algebra is a represen- 
tative of a class of nonlinear algebras which appeared in conformal field theories. It is a product of a quantum version 
of an algebra of area-preserving diffeomorphisms of a 2D cylinder and the abelian Kac- Moody algebra. The presence 
of dynamical symmetry Wi+oo is a general property of gapless one-dimensional systems. The chiral generators Wf^ 
of Wi-i-oo are labeled by the momentum index n,k — 2Tm/ L and by the conformal spin index a = 0, 1, 2, . . . oo. They 
satisfy 

[W:, Wt\ = {pn - am)W:tt^ + q{a, /3, n, ™)W^,?+fr'' + . . . + <5"^<5„+„,ocd(«, n) (23) 

where q{a, /?, n, m) and d{a, n) are known polynomials of their arguments (55| . 

For many interesting cases the Hamiltonian can be written as a linear and bilinear combinations of Cartan generators 
of the Wi+oo algebra |93|. The representation theory of Wi+oo in the Verma module can be constructed in the same 
way as representation theory for the Virasoro algebra (56j starting from the highest weight state and building a towers 
of descendant states. Using this representation theory one can construct the action of Wi+oo generators on bosonic 
Fock space. One can therefore establish a correspondence between Bethe ansatz eigenstates and certain combination 
of states of representation module of the Wi+oo- The evolution operator is factorized in the product of factors 
corresponding to different Cartan generators 

U[t)^e'"' ^\[e'^-'^''' (24) 

where Pa are functions of the interaction strength and a level. When the initial state can be expressed as a linear 
combination of descendant states, the application of U (t) is straightforward. 

In practice it is difficult to deal with the whole nonlinear dispersion, so, for practical purposes one can keep only 
next-to- linear terms in this expansion. In this case the Hamiltonian will include only a finite number of Woo generators. 
For example, in the strong coupling regime, corrections to the Hamiltonian to the first order in the inverse power of 
the interaction strength, 

^ = + ^" ) + + W^] + -.- (25) 

The advantage of this approach is simplicity of including higher-order corrections originating from nonlinearity of 
dispersion, and the potential to solve the problem exactly. The disadvantage of this method comes from its limitation 
to treat gapless systems only. Details of this approach together with several related questions will be presented in a 
separate work [57|. 



E. Form factors and decomposition of the initial state 



This is the most direct approach which we will use further in this paper. In this approach we decompose the initial 
state in terms of the eigenstates of the Hamiltonian. The latter form a complete and orthogonal set of states. To 
compute time evolution of the correlation functions we have to combine several ingredients: (i) find complete basis of 
many body wave functions; (ii) know exact eigenvalues; (iii) determine matrix elements (or form-factors) of various 
operators in the basis of exact eigenfunctions; (iv) find decomposition of the initial state in the complete set of exact 
many-body states; (v) develop effective procedure for summation over intermediate states. 

Ingredients (i)-(ii) are known from the Bethe ansatz exact solution. Matrix elements (iii) were computed in many 
exactly soluble models on the basis of the so-called determinant representation. Recent progress in computation of 
the matrix elements of operators in the basis of the Bethe Ansatz wave functions makes it possible to advance in 
this direction. Decomposition of the initial states (iv) over the complete basis depends significantly on the concrete 
nature of the states. We were able to evaluate these overlaps for a ID Bose gas at large interaction strength for 
two types of initial conditions: all particles in a zero momentum state and all particles in Gaussian wavefunction in 
real space. Currently we perform summation over intermediate states numerically. This imposes certain restrictions 
on the number of particles in our system. We checked that in equilibrium this number is sufficient to saturate any 
correlators to their values in the thermodynamic limit. To illustrate this observation we plot the g2{x), computed in 
the ground state using our method, and compare it to the known analytic expression in Fig.([T]). 

Here we make general comments on the structure of the phase space of form-factors. For a systems with a gap 
only a small number of states needs to be taken into account. Contributions of many-particle states are suppressed 
by the gap. This is the case for the NLS system with attraction, which can be identified with the nonrelativistic 
limit of the quantum sine-Gordon model for which previous statement is also correct (see e.g. review and refs. 
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therein for direct comparison of different contribution to the form-factor expansion of the correlation functions). 
Contributions coming from higher order terms grow in the hmit of long time evolution, but should only generate 
subleading corrections. The situation is very different for gapless systems, such as the NLS system with repulsion. 
Many-particle contributions are suppressed at most as a power law. Hence ideally we need to take into account the 
entire phase space. This makes the problem of summation over intermediate states very difficult. In this paper wc deal 
with the most difficult case of repulsively interacting Bose gas in the regime of strong interaction, when, in principle, 
the allowed phase space is huge and multi-particle states are in general not suppressed. 

The form-factor approach has been successfully applied to computation of equilibrium correlation functions in many 
models 58]. A review of the applications of this method to massive models is given in Ref. [s^. The gap in the 
spectrum leads to the rapid convergence of the form-factor expansion. It is enough to take contributions from a few 
particles to saturate any correlation function. This is not the case for massless models where contributions of multi- 
particle processes are essential. Recently this approach was applied to compute equilibrium correlation functions 
of ID Bose gas in (G^IjUII, (see also Ref. [gl] for the attractive case @|). Earlier studies reviewed in [1] allowed 
to evaluate various asymptotics of the time-dependent correlation functions at equilibrium (for concrete example of 
non-zero temperature see e.g. [63l|.[63|.) 

One of the advantage of the direct application of the form-factor technique is the possibility to consider non- 
equilibrium dynamics of weakly non-integrable systems. In Refs. j65| the form-factor perturbation theory has been 
developed for a class of models which deviate weakly from an integrable " fixed point" . This perturbation theory is 
unusual because its unperturbed states are highly correlated states of interacting integrable theory. We expect that 
this formalism can be extended to treat non-equilibrium dynamics as well. 

The form-factor approach supplemented by a procedure of decomposing the initial state into a complete set of 
many-body eigenstates is universal and can be applied to a large class of integrable models. This method allows 
several extensions and modifications. In the rest of this paper we focus on this method and demonstrate that it is a 
powerful tool for analyzing non-equilibrium evolution of correlation functions. 

In the next section we also show that one can use the intertwining operator to find an overlap coefficient between 
interacting and noninteracting states. We do this explicitly in the strong coupling limit. [9^ 

To conclude this section we note that other possible app roaches to computation of the correlation functions for 



systems out of equilibrium have been discussed in Refs. 66[ and [6 



III. EXACT APPROACH TO NON-EQUILIBRIUM DYNAMICS BASED ON FORM-FACTORS 

In this section we describe the application of the method of form-factors to the computation of non-equilibrium 
correlation functions. This method is general for any integrable system. In this paper we focus on the ID Bose gas. As 
discussed earlier the treatment of gapless system is more challenging because of the large number of contributions from 
intermediate states which need to be included. To overcome this difficulty we devised a special numerical procedure 
of summation over intermediate states. This procedure is described in the next section. 

A. Formulation of the problem 

To fix notations we consider the nonlinear Schrodinger Hamiltonian on a finite interval [0, L] 

H = ( dx[d^-^\x)d^^{x) + c^'^[x)^\x)-^{x)-^{x)] (26) 
Jo 

where the commutation relation between bosonic field "^{x) and its conjugated \['^(a;) as well as the pseudovacuum 
state |0) are defined as Q 

[^{x),^\y)] = 5{x-y), [^{x),^{y)] = [^\x),^Hy)] = Q (27) 
*(a;)|0) =0, (0|*^ = 0, (0|0) = 1. (28) 

The total number of particles and the momentum operators are conserved quantities 



N= / *1'(x)*(x), P = [^\x)d^-^ -{d-^Hx))-^{x)] (29) 

[H, N] = [H, P]=0 (30) 



9 



In each A''-particle sector this system is equivalent to the ID Bose gas with contact interaction potential (Lieb- 
Liniger model) Q, In connection to the previous section we note that inverse scattering transform for this 

problem has been developed in many papers (see [S] , [s!] , [s^ , [i3l for extensive reviews). 

The problem we are asking here can be formulated as follows: suppose we prepared a system in a certain initial 
A''-particle state |V'o^^) (e-g- coherent state can be considered as a superposition of states with different N). This 



state evolves according to the Hamiltonian ([26]). The question is to compute the correlation functions of the following 
type 

{4''^\^\x^,h)^>{x,,h)\4'''>), (31) 
{^//,''^\^Hxuh)^{xuti)^H^2,t2)^{x2,h)\4''^), (32) 
(V'r^|exp(aQ(a;))|4^)), (33) 

where the last formula expresses a generation function for the density-density correlator, 

Qix,t)= f ^\y,t)^iy,t)dy (34) 
Jo 

The density is defined as 

p{x,t) = ^''{x,t)^{x,t) (35) 
Assuming that the initial state is normalized, one can insert a resolution of unity several times, 

into expression for the correlation functions. This suggests to define normalized correlation functions, 

(v;rio(xi,ti)o(x2,i2)iV'r) 

= V V V (^r^l{^i^)({^}^|Q(^i'^i)IMA-)(M^|Q(^2'^2)|{t^}^)(M^|V'r^) , . 

m I ; ({a}|{a}>({m}|{m}>(MI{^}> ■ ^ ^ 

{\}n {m}n {viN 

In the following we concentrate on the density-density correlation function. According to our programme we need 
matrix elements of the operators in the complete basis of states. In this basis time and space dependence are trivial: 
using the Galilei invariance (we implicitly assume periodic boundary conditions) the matrix elements which we need 
are given by the expressions 

mMQ{^,t)\WN) = {e-^^r-Pn _ l)e-(<'^4"')^^Q({^};{A}). (38) 
Here the energy and momentum are given by 

N N 

E^-T^^l ^r-E^.- (39) 

i=i j=i 

In the following it will be convenient to use more compact notations 

4°^ = (^'^^KA}^) (40) 
FQm;M) - (M^|Q(0,0)|{4jv) (41) 
ll{A}|| - ({A}|{A}). (42) 

We can write down an expression for the density-density correlator as 

(4^^(^,0^(0,0)14^)) (43) 
9' A^"\A^°^rFQ{{X}■,{^,}){FQ{{^,}■,M)r 



E E E 



dxidx2^j-^ . ; fr ll{A}||||MIIIIMII 

{A|jv IM/iV jl^lN 

X (e-i(Pr'-^r') _ i)e»tiK")-<')(g«.(P<")-Pr') _ i)e»*^(^i"'-4"') (44) 
where * means complex conjugation. 
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B. Bethe Ansatz Ingredients for ID Bose gas 

As discussed above in Section piEp we need several ingredients: a set of exact eigenstates and eigenenergies of the 
model; an overlap of the initial state with the eigenstates; an explicit expressions for the matrix elements of operators 
in the eigenbasis. 

1. BA states 

The BA states are described by a set of N real numbers A which are given by the solution of the BA equations 
[ll,@,i 

A. - A, 27r , ^ iV + 1 . 



A, + -^2arctan^— ^ = j^l...N (45) 



1=1 



Here the quantum numbers Ij are half-odd integers for N even, and integers for N odd. For eigenfunctions rapidities 
do not coincide, Xj ^ Xk for j ^ k. The whole Fock space is obtained by choosing sets of ordered quantum numbers 
Ij > Ik, j > k meaning that Xj > Xk, j > k. 

From solution of these equations we immediately obtain energies and momenta via Eq.(|39p. 



2. Overlaps between BA states 
The overlap between BA states is given by the Gaudin-Korepin formula ^ Q, [3,[7i|, |zl,[zi,[zl,|6l| 

where 

Sjk{{X}) = Sjk 



N 



L + J2KiX,,Xi) 



1=1 



K{X„Xu) (47) 



2c 

(A,-A;y2 + c2 



^(^.'^^) = 7^ TTY-T^ (48) 



3. Matrix elements 

The matrix elements (form-factors) for the current operator are given by the following determinantal expression 



24l,l25|, 161| 



^^q(Ma^;{a}a.) = n (^^^^) 



where 



[1 

nf=i(Aa - Xj±ic) 



yi±) ^ n.^-i(A^a-A,±ic) ^^^^ 



and where K is given by (pS)) . The matrix inside det() has only real entries. This matrix has size N x N. Real 
numbers {A}, {fj,} are solutions of BA equations P5|) . 
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C. Overlaps of the initial states with Bethe Ansatz basis 

Here we study in details two physically motivated examples of initial states. The first one has its origin in the 
condensate physics whereas the second one describes a Gaussian pulse created in special nonlinear media. Before 
going to these examples we make several general observations. 

1. General remarks 

To compute the overlap with the initial state we consider the following (coordinate) representation for the Bethe 
ansatz states of the NSE 

|*w(Ai, . . . Xn)) = / d^'zxNizi, . . . , zn\Xi. Aw)^'t(zi) . . . *t(^^)|o) (51) 

L V TV! J 

where the function xn is given by 

d d 

XNizi,...,ZN\Xi,---,XN) = const Y\_ (o o hc)det[exp(iA-,-Zfe)] 

JV>j>fe>l ^'^ 
N 



^(-1)''^! exp[i ^ ZnXvn] ]^[A-Pi - X-pk - ice{zj - Zk)]. 



(52) 



It is convenient to define the basic wave function 



|vI/)W = ^vI/t(^^)...vI;t(^^)|0). (53) 
In principle, any quantum state can be expanded as 



(54) 



with 



^|a„|^ = l, \fn{xi, . . . ,Xn)\'^dxi . . .dXn ^ I (55) 

n •' 

Intuitively one expects that the evolution of some special initial states can, under certain assumptions about these 
states and about the Hamiltonian, be relatively simple. This is indeed the case for the ID Bose gas. Using the 
expansion Q and contour integration it was shown in [32] that the state defined as 

l^bNit = O))ord - 0ixi >X2> ...> xn)\^\xi)^\x2) . . . ^\xn)\Q) (56) 

is equal to the state R'^ {xi)R\x2) ■ ■ ■ R\xn)\Q) where R{x) = / ^e"^i?(0 (see Section HEl for definitions). The 
calculation of the evolution of this state is therefore straightforward. 

2. The overlap with a delta function in momentum space. 
We first construct a Fourier transform of the basic state: 

|vE'(gi), . . . *(gw)) = ^=i=[ e^''---■^\x^)dx^], ...[f e^'^-'^-^\xN)dxNm ■ (57) 
\/N\L^ Jo Jo 

At the end of all computations one could take all momenta to be equal qi — q2 — . . . qn to obtain a condensate state 

lV'r^(g))c = [*n<z)]^|0). (58) 
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Let us consider the overlap of this state with the Tonks- Girardeau wave function corresponding to the large c 
limit. There are several reasons to focus on this limit: 1) from our basic interest in application to quantum optics 
with strongly-correlated photons the TG case is the most interesting since effects of fermionizations [tB] should be 
most pronounced in this limit; 2) equations for the overlap are relatively simple and amenable to analytical 
calculations; 3) one could use an expansion in powers of 1/c around TG limit to produce results for the correlation 
functions beyond TG limit. 

In the TG limit the wave function takes the following form [l^, [77|.[78t: 

|*TG> = 4= I dxi... ( dxNXTGi{^^}\{^'^})^^\^l) ■ ■ ■ *^(a;A')|0) (59) 



'N\ Jo 

such that {'^tgI'^tg) = 1 and where 







Here 



VN.L i<j^k<N 
2^ r e Z, N odd; 



L ' \ ?ii G Z + i, iV even. 
The overlap therefore is reduced to the following integral: 

c{i,^^\q)\^Ta)^j^ f dxr f dx2... f dxM J] " ^^fe) det(e""(^"-«)) (62) 

where the factor TV! comes from the overlap with the basic state and denominator comes from the normalization. 
The result for the overlap is given by 

c(^r^(9)|*TG) = (63) 

jV!2^ni<j<fc<jv(Aj--Afc) (_^)f (cos(^))f (sin(:^))f , N even; 

i^a=i(A.-g)ni<.<,<Ar(A. + A, -2g) \ (-z)^(cos(f ))^(-sin(^))^, N odd 

. N 

X (_i)i:£.n.exp(^5]A,)exp(-*7VL|) (64) 

This formula can be simplified further — see Sec. IIV (J2] below. 

Straightforward generalization is a construction of a coherent state in fc-space. One can take a superposition of 
\^^^\q)) states with different iV's for real or complex a: 

oo TV 

l«(9))=e-^l"l^5]^|*W(g)) (65) 

We will not consider this state here. 

3. The overlap with the Gaussian pulse 

In the case of a pulse prepared at time t = we consider a state 

r ^ 

|V'r^)G= \[dx,f{xi,...XN)^\xi)...^>\xNm (66) 

i=l 

and choose a function f(xi...XN) in the form of a product of gaussians, 

N 

f{xu...XN) = Y[eM-x^j2a^] (67) 

i=l 

For simplicity, we take = ct, Vi. The overlap with the TG state is given by 

{^TaHr)G = ^(|)"/^-"exp(-^5:A?)n[Erf(^^-^) + zErfi(^)] (68) 
For all practical purposes the Erf function can be replaced by the gaussian. 
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D. Tonks-Girardeau limit 

From our understanding of the equilibrium structure factor of the ID Bose gas [g^] we expect that the allowed phase 
space of the weakly interacting gas is limited to the vicinity of the dispersion curve of Bogoliubov quasiparticles. On 
the other hand we expect that in the opposite case of the TG limit the allowed phase space should be considerably 
larger. Hence non-equilibrium dynamics should be more involved in this case as well. Here we thus focus on this 
interesting limit in more details. Certain simplifications occur in the equations for the matrix elements, which are 
described below. 



1. Reduction of the form Jactors 
To compute density-density function in the TG limit 



we have to know the following form-factor: 



F^{x, t, {A}, = e^*(^-iJ.)^^ / dx^... dxN ^(-1)[^1+[21 exphzx(Ap - ^JiQ)] 



(70) 



N-1 



exp[-i ^ Xa{Xva - AigJ]- (71) 



A limiting procedure for the form factors in this case has been developed earlier [79|,!8C|,[8l|. We use here the fact 
that the TG limit can be obtained as a double scaling limit of the XX chain. Using 

where we use the boundary conditions e"^^^" = e*-^^*" (—1)^+^. Therefore 

F;(x,i,{A},M) = | jJ-^' Pr,/^^ X n =L''iPoSx,+ '~^^l'~'\ l-S,,)) (73) 

[L e ^ f"^, Ai — ^1, . . . AAr_i — ^jv-l, A ^ ^ L 

and otherwise. 

Note that this form-factor is antisymmetric function with respect to any interchange of momenta for both sets {A} 
and {fj,}. So when they are not ordered in the formula above, we have to assign a factor {-l)[P]^+[PU , where [P] IS a 
permutation index. 

2. Delta function subtraction 

When we calculate density-density correlator we need to subtract a contribution of the particle with itself. 
This contribution comes from commutation relations for Bose operators: 

{X\¥{x,)^\x2)'f{x2)'f{Xi)\,^) = {X\¥{x,)'f{xi)^\x2)'f{x2)W) - {X\^\x,)^ {x,)\,,) 5{x, - X2) (74) 

Considering the TG limit and using (j73p we can calculate the corresponding (5-function contribution 

-{X\^\xi)-^{xi)\v)5{xi^X2). (75) 
For the ground state in the continuous limit this contribution is given by 



1 

2^ 



dqe'" (76) 
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where x = xi — X2 and t = t2 — ti. In general case (|69|) delta function yields an additional term 

-Ipo + j (-l)[^l^+[^l-A;,A,e"i(^"-'"')-"i(^-"^") I ^^e"«-'*^'. (77) 

Now we can combine all components together to calculate the density-density correlation function in the TG limit. 
For the ground state we have 



g^'\Ax, At)=p^ + ^ f f e^*(??-«i) cos[(qi - 92)^] f ^ dqe 



+00 

iq'^t—iqx 



\qi I >TTp J\q2\ <-rrp 

In a general nonequilibrium case correlation function can be written as 



(78) 



+Pq— ^ (^_l'jlP]^ + lP]^> j{^j^^(^gixi{^a~Pb)-itiiEx-E^) _^ ^iX2(Xa-p.b)~it2{Ex-E^)-j 



itE„ 



(79) 



where the first term represents 5 — 5 contribution, the second one corresponds to X ~ v (diagonal part), third and 
fourth terms come from the nondiagonal (\^ \x ^ v) part, and the last term is a delta contribution discussed above. 
Here everywhere A differs from p by one filling number only, and p differs from v by one filling number as well. 

Here, the 5-5 part as well as diagonal parts are symmetric with respect to interchange of momenta since these are 
the products of two antisymmetric functions. On the other hand in off-diagonal terms we inserted the sign factors 
in order to ensure that proper symmetry is preserved: since each F is antisymmetric, the nondiagonal part will be 
antisymmetric as well provided that functions Ax and are symmetric. 

The second term (diagonal part) can be simplified further using the following relation 

N N 

E f{K)=J2f{Xa)-J2f(f'>') (80) 

Aqt^/ai .../Ajv a— 1 6—1 

for arbitrary function /. Then 

L L 

^ N N 

+ J2 E E |AA|^e'^-^»-"^^e^"^''+"^^ (82) 

a=l b=l 

This identity can be also used to simplify the nondiagonal part. 



3. Deviation from the TG limit 



The main difficulty of our approach is computing overlaps with the initial state in a compact, analytical form 
for arbitrary interaction strength. At present we do not have a complete solution of this problem for an arbitrary 
state. For particular types of states the problem can be attacked using one of the formalisms of Section II. Thus, 
in particular, if the initial state can be expressed in second-quantized notations via creation-annihilation operators 
acting in momentum space, then the second-quantized version of the intertwining operator can be used. If the initial 
state has a simple form in coordinate space, the coordinate version of the intertwining operator might be useful. One 
form of this operator in real space, mostly useful for the 1/c expansion around TG limit, has the form of differential 
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operator acting on the TG wavefunction [ll,!!^. It relates the TG wavefunction (which is almost fermionic) with the 
finite-c Bethe state 



\BA), 



' n 



d 



l<i<j<N 



c dxi dxi 



|det(e'^"^")|. 



(83) 



Manipulations using intertwining operators will be discussed in more details in the future in connection to concrete 
physical problems. Here we note that expanding the intertwining operator up to order 1/c one observes (96j that a 
correction to the overlap coefficients up to the order 1/c (neglecting the boundary effects) for the condensate initial 
state have the following form 



A 



(0) 



Af(l 



(84) 



Note also that deviations from the TG limit includes other ingredients. One of them is an expansion of rapidities in 
powers of 1/c, 



A, -(2.7 -TV- 1) 



L 



1-2^+4(^)2 



c 



4^2 

T 



;-)' + —N{2f + {N+ i){N - mi^f 



(85) 



Using this expression one can get a systematic expansion of the form-factors. We do not need these 1/c corrections 
since our procedure uses numerical expressions of the form-factors. 

We finally note that a computation of the overlap of the arbitrary BA state (0|C(Ai) . . . C(A„) with the initial state 
of a special form vE'(a;i) . . . \E'(a;„)|0) can be considered as a particular case of evaluation of the multiple point form 
factor if we regard the pseudo vacuum state as the BA state. The problem of computation of such a multi point 
formfactors can be solved using the formalism of multi-site generalization of the NLS problem. This problem was 
addressed in 67|. Another way to go beyond the TG limit is to use different formalisms described in Sec. II. 



IV. NUMERICAL TREATMENT 



A. Introduction 



In this part we present results of numerical computations of correlation functions for various non-equilibrium initial 
states and their time evolution. In order to do that, we needed to sum expression (j44p over all initial A, intermediate 
jjL, and final v states. 

Although we are interested in thermodynamic limit L — oo, — > oo, numerically one can compute correlation 
functions for finite size systems only. In order to infer thermodynamic limit from our finite size results we can set 
the size of the problem [N and L) to be sufficiently large so that numerically computed correlation functions would 
be sufficiently close to the ones for an infinite system. One of the criteria for such sufficiency could be comparison 
of the numeric results to analytic results for the system for which thermodynamic solution for correlation functions 
is exactly known. For example that we already know the exact expression for a correlation function in the ground 
state (j78p . So we can use it to check the accuracy of our numeric procedures. 

Let's estimate the cost of a brute force approach towards summation (|44p for a reasonably large system size. Because 
we are interested in non-equilibrium states, states in A, /x, and v are filled way above the "Fermi momentum" . We can 
estimate the number of terms in the sum (H^ by taking N ^ 100 and momentum cutoff ~ lOA^. Then the number 

of terms in the sum ^ (^^) - 10'*20. Obviously, the sum can't be taken in a straightforward fashion. 

Now we would like to look for possible simplifications of the problem. The first hint comes from [63| , where authors 
claim that when one considers summation of the form 

Y.\{GS\F\n)\\ (86) 

where GS is a ground state, then the sum has a very limited number of major contributions, most of which are one 
particle. Note that this statement is not directly applicable in our case because we are not working with a ground 
state and also because we have a double summation, which might have different major contributions due to larger 
phase space, or not to have such major contributions at all. 

Second simplification might come from the fact that we are interested in a large interaction constant c case. For our 
purposes we can consider the case of large c. As we found in the previous section, relevant quantities become much 
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simpler. The most prominent feature of the TG limit is that the only non-zero matrix elements (A|i^|/x) come from 
zero and one particle processes (j73|) . It means that in this limit for a given A we need to do a two particle summation 
in (|69|) with one particle transition A — ?► /i and one particle transition ^ i/. This observation also tells us that for 
large, but finite c, processes with small number of transitions should dominate. 

Finally, structure of overlap factors A\ for different initial states might provide some clues on how to simplify our 
task. In particular for the initial states we are going to investigate, overlap factors can be simplified analytically so 
that we do not need to perform the summation (j52p over all permutations. 



B. Implementation 

Now we would like to discuss ideas behind concrete implementation of our calculations. Almost everywhere we 
consider the limit of large c, though we can do numeric calculations for arbitrary c at a cost of substantially higher 
cpu time. 

For a ground state we do several large-c calculations to check what values of c can be treated as "infinite" and to 
investigate multiparticle contributions. For large, but finite c, it is possible to further simplify expressions (|46l) - (|50|) 
by doing 1/c expansion, as explained in the previous Section. We implemented both direct and simplified versions of 
matrix elements and overlap calculations and did not find any visible discrepancy in the results, though performance 
improvement was substantial. 

After we sum up the expression 

J2 AxA,{X\F\f,){t,\F\,,) (87) 

over intermediate states /i and final states v for a given initial state A, we are left with a problem of summation over A. 
But the phase space is huge, therefore we can't use a direct summation. We use Monte-Carlo summation instead [83| . 

Namely, we sample set of states {A'} from the set of all possible initial states {A} with probability of state A' to be 
selected being P\. Then we can replace summation in (|87p over A with summation over A' to get 



rMC 



Ex',,,.AyM^'\m{f^\FW)/Px' 



As we increase the size of sample {A'}, the value of this sum converges to the true value of (|87| for any probability 
distribution P\. Depending on the pick of the probability distribution function convergence time can be very different. 
For instance if we take a uniform distribution over the entire set {A}, than we can see from Eq. (|68|) that weights A\ 
can be very different. As a result, many contributions to the (1881) will be negligibly small. On the other hand if we 
are picking A such that \A\ \ is around its maximum, all the contribution to the sum will be substantial. 

Indeed, it is known [83], that the sum optimally converges to the true value if the distribution function is proportional 
to the expression we sum. In our case it would (roughly speaking) mean P\ oc A\ A^, where v differs from A at 
most by two particle process. As a reasonable approximation, we use Pa (x A\. 

The next question is how to sample such states A'. We cannot sample states from the distribution directly, but we 
can use Gibbs sampling [s^ to do that approximately. Note that any overlap factor can be expressed as a function of 
state A and parameters of the problem. A itself is a function of a vector of quantum numbers / (|^5t . We can fix all the 
components of / except one {Ij). By varying this component and resolving the equation (1451) relative to A, we can find 
conditional probabilities P{Ij\Ii . . . Ij-ilj+i . . . /„) oc , where is a solution for A given quantum number Ij 
(and the rest of quantum numbers, which stay the same). This way we can sample Ij from the P\ distribution given 
all other quantum numbers li^j are fixed. Then we iterate the procedure with index j running from 1 to A'^ several 
times. Using this algorithm we obtain a set of statistically independent states pooled with probability oc A\ (83 |. 

In theory one can do Monte-Carlo summation not only over A states, but over ^ and v as well. Nevertheless, 
we found such approach impractical because of a much slower convergence stemmed from very poor cancellation of 
various harmonics. 



C. Results 



1. Ground state 



First we would like to focus on correlation functions for the ground state. We have several reasons for that. First, 
exact analytical solution for such correlation functions is known (|78p , hence we can at least partially verify the validity 
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FIG. 1: Upper panel: calculated and theoretical spatial correlation function for a ground state for c = 100, c — 1000, and 
c — >■ oo. Lower panel: zoom of the upper panel. 



of our approach. Second, we do not have to sum over A and v, so computations are very fast and precise. Therefore 
we can use this case to investigate the effects of our approximations — what finite values of interaction constant c 
can be considered as infinite, how one particle approximation affects results, and what are the finite size effects of our 
calculations. 

First we looked at spatial correlator, calculated at L = 51, N — 51. For the ground state it is translationally 
invariant both in space and time. On Figs. [l]we show a theoretical correlator along with ones calculated at c = 100, 
c = 1000, and c = oo for one particle processes only. We see that for c — 1000 the correlator is almost indistinguishable 
from the theoretical c — >■ oo limit, though c = 100 slightly deviates from that limit. Hence all our further c — >■ oo 
results apply to c > 1000 case as well. Also note that theoretical correlation function perfectly overlaps with one 
calculated for c = oo, hence our choice of system size can be considered as a thermodynamic limit. 
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FIG. 2: Theoretical and calculated temporal correlator for a ground state in the c — oo limit. Real, imaginary parts, and 
absolute values are shown. Only three curves are visible because theoretical and calculated plots perfectly overlap. 

As a separate run we calculated the same correlators for the same finite values of c while taking into account all 
two-particle contributions and major three-particle contributions in addition to the original one-particle ones. We 
found that these contributions are negligible. We conclude that it is safe to ignore them in the large c limit. 

On Fig. [5] we show a temporal correlation function for a ground state for c — > oo as well as a theoretical value given 
by (j78l) . We see that they perfectly overlap except for the tiny region around x=0 because of momentum cutoff in 
our numeric calculations. 

2. Delta function in momentum space 

We consider the case when before interactions are switched on all particles are initially in the same state with a 
given momentum p. Such state can be prepared experimentally from the condensate at rest using a Bragg pulse. 

Before we proceed with calculations, we would like to take a closer look at overlap factors A\ for this state ([M]) . 
We are starting with the equation for the projection of our state onto eigenfunctions (TG) of the Hamiltonian 
{'^'^^\q)\^TG)- Ignoring coefficients independent of the TG state and non-singular at g = 0, we have for odd [97|] N 



vI/W(g)|M'TG)(x 



<j<k<NV^] 



{\j - Afc) 



ntr(A. 



sm 



Lq 
2 



(89) 



Note that our system is discretely translationally invariant in k space: we can shift all and q by the same constant 
for arbitrary n, change quantum numbers Ij by n, and all the equations (|45p - (|64p will hold. 

Neglecting finite size effects (resulting in discretization), we can choose n such that q ~ Then we shift all A.^ 

by q, and consider the overlap ((89|) factor in the g — ?> limit. 

Sinus in the expression (j89|) results in a ^^'^ order zero in the numerator in g — > cx) limit. Zeroes of the denominator 
at q — are given by one of Xi = in the first product, and all the pairs < j such that A.^ = — Aj in the 

second product. Maximum order of zero we can achieve in the denominator is when we have one zero Xi, and for each 
of the rest of lambdas we have one exactly opposite to it. For such configuration zero's order of the denominator is 
also . In this case they cancel resulting in a non-zero contribution. Any other configurations will result in the 
lower order of zero of the denominator, hence zero overlap. We conclude that all the states with non-zero overlap are 
symmetric. One can easily repeat this calculation for even TV and show that states in this case are symmetric as well. 

We can draw an immediate conclusion about the structure of the sum (|88p we use to calculate correlation function. 
Remember that for an infinite c {X\F\fi) has non-zero elements only if states A and /i differ by at most one quantum 
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FIG. 3: Three-dimensional image (coordinate, time and intensity) of the correlation function for the initial state corresponding 
to the delta peak in momentum space. 



number (|73p . Therefore in the summation ()88|) A is symmetric, fi has one quantum number changed relative to A, 
and ly, being symmetric, has one quantum number changed relative to /i. But it is possible only if either |A) = 
or if the changed quantum number in |A) is symmetric to the corresponding change for — \i>). I.e. if first 

we change li /{, than later we should change either — )■ or —li = Ij — )■ /j = — As a result for a given |A) a 
selection of allows only two possible states for {v), strongly reducing the phase space we need to consider. 

Using the symmetry we can also simplify the expression for overlap factors (j64p . Because in the sum (j88p normaliza- 
tion of A\ does not matter, we drop all the constant terms and can consider ((89)) instead. Let's denote n — {N — l)/2 
(N is odd), and indexes run from —n to n. Then we can rewrite the numerator as 

n {Xj - XkfiXj + Xkf n -2^' (90) 

0<j<fc<n 0<j<n 

and the denominator as 

(-g)^2^ n (A,-Afc)2(A,+A,)2 H (91) 

a<j<k<n a<j<n 

so that (|M)) simplifies to 

Ax « \ ■ (92) 
11aj>o -^j 

In case of interest one can easily restore the constant coefficient and also prove that this expression is valid for even 
N as well. 

Using the conclusions drawn above about relative structure of A, /i, and v states, and using expression (|92p for 
overlap factors, we perform a Monte Carlo summation over states for N = 101, L = 101. On Fig. |4]we show a spatial 
correlation {p{xo, t())p{xo + x, to)) as a function of x for xq = and different values of Iq. Basically the plot represents 
time evolution of x correlator starting from the moment we momentarily switch on a very strong interaction. We 
can see that at the moment to = the entire correlation function is not distorted except for x = 0. This is because 
when we immediately switch on the interaction, particles stay where they are. They just don't have any time to move 
anywhere, hence the original non-interacting correlation function, which is flat, is preserved (small distortions around 
a; = are artifacts of our calculations; if to take into account their errors, they are indistinguishable from 1). But as 
the time goes by, repelling interaction generates a wave, which transfers the matter away from the particle sitting at 
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FIG. 4: Left: spatial correlator for the initial state corresponding to the delta function in momentum space for different to. 
Right: zoom of the left panel. 
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FIG. 5: Real and imaginary parts of a temporal correlator for the initial state corresponding to the delta function in momentum 
space. 



X = 0. On the figure we observe the evolution of this wave. When the time becomes sufficiently large (i > 1), the 
correlation function stabilizes and does not evolve anymore. 

We also calculate a temporal correlator (/5(a;o, to)p(a^o, + t))- Figure [H] shows real and imaginary parts of the 
correlator. Because the state is translationally invariant, correlator does not depend on xq. Interestingly, it also 
almost does not depend on tg. All the plots for different overlapped, and the difference was invisible. On this plot 
we also observe two "effects" , which are artifacts of our calculation scheme because of introduced momentum cutoffs: 
oscillations in the imaginary part and real part not going exactly to zero at t = 0. 



3. Gaussian pulse 



Now we proceed with the initial state — a Gaussian pulse in space, and hence in momentum space. Overlap factors 
for such state are given by the expression We neglect finite size effects by considering the limit L ^ a, for which 
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FIG. 6: Evolution of the correlation function for the initial state corresponding to the Gaussian in the real space. 



expression ()68p . up to a constant factor, can be written as 

^ACxexp(^-^^A,?j. (93) 

In this case we do not have a nice symmetry we had in the previous section, so we have to sum over ah one particle 
/i states and two particle v states using expression (j88p . Summation over all two particle processes substantially 
increases the computational time, so in order to partially compensate for it we used smaller values for problem 
parameters i = 21 and N = 21. For a packet width parameter we use value a = 15/(27r). 

On Fig. [7]we present time evolution of spatial correlator {p{xQ,to)p{xo + x, to)) for xq = and different fixed times 
to. We also show errors of our calculations to distinguish true correlation function features from noise. 

We would like to emphasize several features visible on these plots. First, for to = 0, as expected, we see a dip 
at X — 0. Interestingly, the width of this dip is finite and is not determined by the momentum cutoff. From 
computational point of view we can explain it by the fact that the overlap factor ([M)) for Gaussian state strongly 
inhibits high harmonics, therefore we do not have momenta high enough to make this dip infinitesimally narrow. This 
is a finite size effect because of a small number of particles. As parameters L, N, and a increase proportionally, width 
of the dip goes to zero. Second, oscillations of the correlator at to = around the dip xq = are not determined by 
the cutoff and much bigger than the error of our calculations. This is again the effect of a finite size system. Third, 
as we've already seen in previous section, for times ^ 1 the correlation function stabilizes. But the process of 
stabilization is quite remarkable. We do not have a widening of a Gaussian peak. Instead it looks like the peak has 
the same width, but its height decreases with time. Probably, qualitative explanation lies in the fact that high density 
regions have high interaction energies, hence in these regions we have fast particles which quickly leave the region. 
I.e. we have a ballistic scenario, which differs from the diffusion scenario for which we would expect a widening of a 
Gaussian peak with time. 
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FIG. 7: Spatial correlator for different to for a Gaussian initial state and its error. Two lines for each value of to lie one standard 
deviation above and below relative to the average value. 

Unfortunately, due to some peculiarities of computations we had to conduct, we were unable to produce reliable 
results for a temporal correlator. Results, though reproducing an overall shape of the correlator, were very noisy 
because of beats attributed to non-exact cancellation of various harmonics. We do not provide a plot for these results. 

V. DISCUSSION AND SUMMARY 

In this paper we addressed a problem of non-equilibrium time evolution of one dimensional Bose gas with contact 
interaction from the general perspective of dynamics of integrable systems. Several approaches have been proposed 
recently for analyzing time evolution of integrable many-body models, including Quantum Inverse Scattering method, 
the formalism of Intertwining Operator, and Extended Conformal Symmetry. After critically reviewing these ap- 
proaches we concentrated on conceptually the simplest method, based on using Bethe ansatz to decompose initial 
states into precise eigenstates of the interacting Hamiltonian and using form factor expansion to calculate time evo- 
lution of correlation functions. The main difficulty of this method is that it requires summation over a large number 
of intermediate states. In this paper we focused on the regime of strong repulsive interactions between bosons and 
developed an efficient numerical procedure for performing summation over intermediate states. We analyzed two 
types of initial states: all particles having zero momentum and all particles in a gaussian wavepacket in real space. In 
both cases we find a non-trivial time evolution of the second order coherence g2{xi, X2,t) = {p(xi,t), p{x2,t)), which 
reflects the intrinsic dichotomy between initial states and the Hamiltonian. Initial condensate states exhibit bunching 
at short distances, whereas strong repulsion in the Hamiltonian introduces strong antibunching. 

For the initial state which has all particles in a state with zero momentum, results for 52 are summarized in Fig. |31 
Antibunching at shortest distances is present at all times reflecting repulsion between particles. A more striking 
feature is the appearance of bunching and oscillations in 52(2^1 — X2,t) at intermediate time and length-scales. At 
transient times 172 has Friedel like oscillations as a function of X1—X2, which disappear at longer times. Crystallization 
of the TG gas in the process of non-equilibrium time evolution was also discussed in Rcf 20] . At longer times the form 
of g2 is qualitatively similar to what one finds for the equilibrium Lieb Liniger model with intermediate interactions 
at high temperatures [s^ . Fig. |3] also provides some support to the idea of light cone formation discussed in Ref. [s^ . 
Some of these features are similar to quench dynamics in other integrable systems [13, [s^ . 

For the initial state which has all particles in a gaussian wavepacket in real space, time evolution of (72 is shown in 
Figs. 16171 As in the previous case we observe antibunching at the shortest distances. For intermediate time scales it 
is followed by a pronounced peak in 172 refiecting dynamic bunching. This system would be naturally characterized 
by the time dependent short range correlation length, which increases with time, reflecting expansion of the system in 
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real space. Oscillations in (72 which we observe in this case are small and are most likely related to the finite number 
of particles used in the analysis. Hence transient states of this system would be more natural characterized as a liquid 
rather than a crystal. 

Predictions made in our paper for the behavior of g2 should be possible to test in experiments with ultracold atoms 
and photons in strongly non-linear medium. 
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VII. APPENDIX A: INVERSE SCATTERING TRANSFORM AND ALGEBRAIC BETHE ANSATZ - 

BASIC CONCEPTS 



The nonlinear Schrodinger Hamiltonian represents the simplest example of a system solvable by the algebraic Bethe 
ansatz [ss", "gs"! . For review and many details see 4] . Here we overview the basic construction with connection to the 
inverse scattering transform in order to fix notations and for the sake of self-consistency. 

The Zakharov-Shabat method starts from transforming the field '^{x,t — 0) of nonlinear Hamiltonian ([T]) at time 
t = into a set of " scattering data" given by the following linear problem 

i^<^{x,i) ^ Q{x,0^{x,0 (94) 



where 



Q(x,0- ( , "^f^"^^ ) (95) 



X 



The solution <I>(a;, ^) of this equation is defined by the condition |\I>(2;)| — > as x — >■ ±00 and by the properties of Jost 
solutions, 



(97) 



Rewriting the linear equation above with the boundary conditions as an integral equation gives the Gelfand-Levitan- 
Marchenko equation 



ei(x,Oe'^«/' = -N/S / dy0(2/-.T)e^€^^'t(y)x2(2/,Oe-'^«/', (98) 

•/ — 00 

6(x,0e-"«/'-l + *\/^ / rfy0(2/-x)e-'«^xi(2/,O*(y)e'"*/' (99) 



which can be solved iteratively in powers of ^/c 

00 n „ n 

Xi(x, Oe"«/' = ^ c" II dx^ll dy,0{xo > yi > • • • y„ > a;„ > x)e''«(^-o -^'^^1 v.) (lOO) 

n=0 i=0 j = l 

X ^Hxo)---^\xr,)^{yr,)---^{yi), (101) 

00 n— 1 ^ n 

X2{x, Oe-"«/2 = 1 + E ^" n dx^ll dy,e{xo > yi > ■ ■ ■ x^-i > 2/„ > x)e^«(^-o' "--^"-i (102) 

n=0 4=0 j = l 

X 1't(a;o)---*^(x„-i)^'(y„)---^'(yi) (103) 

where 9{xi > X2 > ■ ■ ■ > x„) — 9{xi — X2)0{x2 — x^) ■ ■ ■ 9{xn-i — Xn)- This solution can be written as a solution for 
the scattering data, A{£^) and B{^). Finally, defining the refiection operator as 

R{0=mO]-'B{0 (104) 
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one obtains the expansion ([2]). Operators R{^) satisfy the Zamolodchikov-Faddeev algebra ([6]). In this formahsm, the 
Bethe ansatz states are constructed by appUcation of Zamolodchikov-Faddeev operators to the pseudovacuum state 
|0), 

\^N)BA^R\Ci).--RHUm (105) 

These states are complete for the case of repulsion. Using the expansion ^ one can show |32|] that the special ordered 
initial state of the form 

\Mord = 0{XI >X2> ...> Xn)\-^\xi)¥{x2) • . • *^(xjv)|0) (106) 

is equal to the following state 

R\xi)R^X2)...RHxn)\0) (107) 

where 

Rix) = J §^^"^m- (108) 

We note that the evolution of this state can be obtained explicitly, as discussed in Section II. The case of attraction 
is discussed in Appendix C. 

The transfer matrix is a central ingredient of the construction of the inverse scattering transform on the finite 
interval on the lattice (one performs a space discretization procedure for the LL model first). It can be represented 
as a matrix 

The transfer matrix is constructed as a product of monodromy matrices Lj{^) for each site j, T{^) = Ylj Lj{£,), which 
in the case of NS problem has the form of matrix Q which is one of matrix forming the Lax pair. Trace of the transfer 
matrix commute for different ^ and is a generator of integrals of motion. The integrability condition is expressed by 
the special property 

R{( - emo ® TiC)) = ine) ® mmi^ - e) im 

where the matrix _R(^, ^') is a solution of the Yang-Baxter equation, which for the case of NS problem belongs to the 
class of 6-vertex model 

i?12(e)^13(e + /i)i?23(M) = i?23(M)^13(^ + /i)i?12(0 (HI) 

Rab = f3Iab + aPab, a ^ ^ ^\ ^\ , 13 ^ ^ ■ (112) 

t,a- t,b- iC t,a- t,b- 1C 

where lab and Fob are identity and permutation operators acting on the tensor product of two single-particle spaces 
indexed by a and h. In the finite size lattice formalism the diagonal entries of the matrix (more precisely its trace 
t(^) = TrT(^) = A(^) -I- £'(^)) generate a conserved quantities of the model whereas the off-diagonal elements i3(^) 
and C(^) act as creation and annihilation operators of pseudoparticles. The n-particle Bethe ansatz eigenvectors are 
constructed using these off-diagonal elements of the transfer matrix as 

|*„({A„})) = B(Ai)B(A2) . . . B(A„)|0) (113) 

whereas the bra- vectors are constructed as 

(vI/„({A„})| = (0|C(Ai)C(A2) . . . C(A„) (114) 

The algebraic Bethe ansatz deals with diagonalization of the trace r(^) of the transfer matrix, thus solving the 
eigenvalue equation t[X)\BA) = A(A, {Aa})|i?A). This can be done using the commutation relations between operators 
A(A), i?(A), C(A), I?(A) which can be deduced from the Yang-Baxter algebra. This procedure leads to expression for 
the A(A, {Ao}) = a(A) 0"=! /('^ ~ ^a) + /?(A) Ha^i /('^a ~ •^) which can be further related to the energy and to the 
consistency relation, known as a Bethe ansatz equations for the momenta Aa of the pseudo particles, 

Oii)^a) 77 Afc - Aa 
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where the function /(A) is determined by the elements of the i?-matrix. The connection with the continuum version 
of the NLS model is established by a limiting procedure of lattice size going to zero while keeping the density constant 
on a finite interval. 

Due to the work of Sklyanin [88| it is possible to formulate even more general NS problem which includes the 
boundaries. The boundary problems can be divided into two categories, soliton-preserving and soliton-non- preserving. 
For recent progress in solution of the boundary NS model see [89| . 
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The graph F with q vertices and p oriented edges is called collision graph if a) for every vertex there is an edge coming 
into it or going out of it; b) there is at most one edge between any two vertices; c) all vertices can be labeled in an ordered 
way from 1 to q such that edges go from smaller to larger number. 

It is also interesting to note that later on the same authors constructed similar second-quantized approach to the fermionic 
Yang-Gaudin model as well 

The authors are grateful to M. Zvonarev for discussions on this issue. 

Cartan generators of any algebra are those which commute with each other. In our case it means that [Wq , Wq] = 0, Va, /3 
Note that direct quantum-mechanical computations has been done in Ref. [l3|; they agree with up to multiplicative 
factor. 

[95] We also comment on some particular type of models where these and other overlaps can be found explicitly for various 
coupling constants. This is the case for some type of Gaudin models, e.g. Richardson model of mesoscopic BCS pairing and 
the central spin problem. In these models the BA states do not depend explicitly on the strength of the coupling constant, 
so one could use the Gaudin-Slavnov determinant formula for these overlaps. 

[96] We note that n^.i Eti det(e-"-)d.^ ... d.i fefn!^ -l^^i^M 

which is equal then to iNJ2'^^^hAf\ where Al = /o°° /q"' ■■• /o"""' 11^=1 e"'""" det(e*""")da;jv ... da;i '=° 

ni<i<j<w(fci-fcj) 
n^=i k„ni<i<j<N(ki+kj)- 
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[97] These calculations can be easliy repeated for even N resulting in the same conclusion. 
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In this paper we study nonequilibrium dynamics of one dimensional Bose gas from the general 
perspective of dynamics of integrable systems. After outlining and critically reviewing methods 
based on inverse scattering transform, intertwining operators, q-deformed objects, and extended 
dynamical conformal symmetry, we focus on the form-factor based approach. Motivated by possible 
applications in nonlinear quantum optics and experiments with ultracold atoms, we concentrate 
on the regime of strong repulsive interactions. We consider dynamical evolution starting from 
two initial states: a condensate of particles in a state with zero momentum and a condensate of 
particles in a gaussian wavepacket in real space. Combining the form-factor approach with the 
method of intertwining operator we develop a numerical procedure which allows explicit summation 
over intermediate states and analysis of the time evolution of non-local density-density correlation 
functions. In both cases we observe a tendency toward formation of crystal-like correlations at 
intermediate time scales. 

I. INTRODUCTION 

Quantum many-body systems in low dimensions can not typically be described using mean-field approaches. This 
makes analysis of their non-equilibrium dynamics particularly challenging. However following demand from recent 
experiments, certain progress has been achieved in developing theoretical methods which can address this problem. In 
this paper we use Bethe ansatz solution to study nonequilibrium dynamics of the one dimensional Bose gas interacting 
via contact interaction. This microscopic model of i5-interacting bosons is called the Lieb-Liniger (LL) model 
It belongs to a general class of models which can be studied using Bethe ansatz solution [1] . These exactly solvable 
models are characterized by an infinite number (in the thermodynamic limit) of conserved quantities which originate 
from the nature of collision processes. In the LL model the two particle collisions can not change momenta of scattering 
particles but only give rise to a phase shift. Moreover any many-body collision is factorizable into a sequence of two- 
body scattering events. While implications of the exact solution for thermodynamic properties have been discussed 
before [U, HQ, their consequences for non-equilibrium dynamics are not known. This will be the central question of 
our paper. 

Systems of one dimensional Bose gases with contact interactions have been recently realized using ultracold atoms 
[Toj and current experiments allow wide control over parameters of the microscopic Hamiltonian (see [ll[ , for recent 
review). For example, the effective strength of the repulsive interaction can be tuned either by changing the density 
of the atomic cloud, or by modifying the strength of the transverse confinement, or by changing the scattering length 
using magnetically tuned Feshbach resonance. In particular a very interesting regime can be achieved when repulsion is 
so strong that bosonic atoms become essentially impenetrable 13 1 and the system undergoes effective fermionizationf^ 



|8| . Recent studies of thermodynamic properties demonstrated excellent agreement between experiments and the exact 
solution [1, . However parameters of the Hamiltonian can also be changed dynamically. In this case one needs to 
study dynamics starting with the initial state which is not an eigenstate of the Hamiltonian. Generally such dynamics 
• y—{ involves superposition of coherent evolution of all eigenstates of the system. This is the problem that we address in 
this paper. 

^ Another class of systems where the LL model appears naturally is light propagation in one dimensional photonic 

- - ■ fibers. It has been known since the sixties (see e.g. ref [l3|) that propagation of classical pulses of electromagnetic 
waves in one dimensional nonlinear medium can be described by the nonlinear Schroedinger equation, which can be 
interpreted as the operator equation of motion arising from the LL microscopic Hamiltonian. In this case nonlinearity 
is proportional to the nonlinear polarizability of the medium, x^*^^ (see e.g. Ref [l5?]). Typical propagation problem in 
one dimensional nonlinear fiber corresponds to the classical limit of the quantum LL model subject to nonequilibrium 
boundary and/or initial conditions [16]. Earlier analysis assumed that photon nonlincarities come from interaction of 
photons with two level atoms, in which case nonlinearity corresponds to the attractive interaction in the LL model [l7| . 
Corresponding quantum models have bound states which eventually lead to formation of solitons and classical regime 
[l^ . Earlier analysis of optical systems also assumed regime of weak nonlincarities since increasing nonlincarities 
in two level systems would lead to strong photon losses. However recent techniques utilizing electromagnetically 
induced transparency [T^ make it possible not only to realize effective repulsion between photons but also to increase 
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dramatically the strength of nonlinearities(20j. Optical realizations of the LL model should allow direct measurements 
of both the first and the second order coherences [21[, g^^\t, r) = {"^'^t)"^!{t + T)) and g^'^\t, r) — {p(t) p{t + t)) , where 
'i'{t) is the amplitude of the propagating electromagnetic field and p{t) — is the field intensity, measuring 

the number of photons. 

There are three main reasons why we undertake a detailed analysis of the non-equilibrium quantum dynamics of 
the LL model. The first reason is that dynamics of integrable models, such as the LL model, should exhibit special 
features originating from an infinite number of integrals of motion. One manifestation of constraints on the phase 
space available for dynamics should be the absence of thermalization and ergodicity ^ ,22,. .23,] . Conservation laws 
originating from integrability also lead to the so-called Mazur inequalities [l^l , which should result in special transport 
properties. For example, the possibility of an infinite Drude weight has been discussed for several integrable models 
of lattice fermions [25|. Our second motivation for analysis of the LL model is to use it as a testing ground for 
developing new theoretical methods and techniques, which can then be applied to a broader class of models and 
systems. Our work should provide a general framework for understanding dynamics of integrable systems and for 
future developments of non-perturbative methods in the study of non-equilibrium dynamics. We emphasize that our 
analysis uses special properties of exactly solvable models, and the methods we apply here are fundamentally different 
from the more conventional perturbative approaches discussed in Refs. (26j.[27j. Finally we point out that recent 
progress in the areas of ultracold atoms, quantum optics, and low-dimensional strongly-correlated materials makes it 
possible to fabricate concrete physical systems, which can be accurately described by the LL model. Thus our third 
reason for analyzing this model is that theoretical predictions for the time-dependent evolution of correlation functions 
can be measured directly in experiments. Specific system in which one can realize dynamical experiments discussed in 
our paper are one dimensional Bose gases in optical lattices and magnetic microtraps (see Refs. [l^H^ for a review). 
Realization of such experiments should provide a unique opportunity to study non-equilibrium dynamics of strongly 
correlated exactly solvable systems. 

In this paper we calculate time evolution of the density-density correlation function for two different condensate- 
like initial states: when all particles are in a state with zero momentum and when all particles are in a Gaussian 
wave packet in real space. There are two main reasons why we focus on these initial states. Firstly, both of these 
states correspond to important physical systems. A condensate of particles in the zero-momentum state provides a 
natural initial state for the discussion of dynamics of a Bose condensate when interaction is suddenly switched on. 
The Gaussian wave packet in real space is a prototypical example of the photonic pulse in quantum optics. Secondly, 
considerable progress can be made in the analytic calculations of the time evolution arising from these two states. 
Also numeical calculations can be done quite efficiently. We emphasize here that it is not the case for a generic initial 
state. Understanding dynamics of the generic initial state requires developing additional theoretical tools. 

This paper is organized as follows. Section II provides a brief critical overview of several different methods which 
can be used to describe non-equilibrium dynamics of integrable quantum systems such as the LL model. Section 
III provides an in-depth discussion of one of these approaches, the so-called form-factor technique, which relies on 
numerical summation over intermediate states. The overlap with the initial states in the strongly interacting regime 
is found using the method of intertwining operator. In section IV we focus on the non-equilibrium dynamics of the 
many-body system described by the LL model. There we choose two particular initial states as explained above. 
Summary of our results and conclusions are given in Section V. To provide additional background to the readers, in 
appendices we provide basic facts about the algebraic Bethe Ansatz formalism and the inverse scattering transform 
for the model of one dimensional Bose gas with both repulsive and attractive interactions. 

II. A REVIEW OF METHODS TO INVESTIGATE NONEQUILIBRIUM DYNAMICS OF 

INTEGRABLE MODELS 

In this paper we analyze the time evolution of correlation functions in integrable models when the initial state is 
not an eigenstate of the Hamiltonian. Here we give an overview of (some of) the possible approaches to studying 
nonequilibrium dynamics of integrable systems. We only discuss approaches which explicitely use the property of 
integrability and exclude more conventional techniques which rely on approximations and perturbative expansions. 



Here ^'(a:) is a bosonic field in the second quantized notations, c is an interaction constant, L is the size of the system. 
The equation of motion for this Hamiltonian is the Nonlinear Schrodinger (NS) equation. The first quantized version 





(1) 
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of this Hamiltonian corresponds to the Lieb-Liniger model of Bose gas interacting with the contact interaction. The 
Bethe ansatz solution of this model has been given in refs. Q , [1] ) Hi > dl • 

Below we discuss several possible approaches to the study of nonequilibrium dynamics of nonlinear integrable 
models, such as the LL model ([T|). 



A. Quantum inverse scattering method 

One approach to analyzing dynamics of model ([1]) would be to use the formalism of the inverse scattering problem. 
This method relies on the solution of the quantum version of the Gelfand-Levitan-Marchenko equation (see Appendix 
A). From the Inverse Scattering Transform (see e.g. [3ll | for review) it is known that the solution is given by the 
following infinite series representation 

OD N N 

*(^) = E / n $ n ^9Ni{p}, {k}; x)rHpi) ■ ■ ■ RHpN)RikN) ■ ■ ■ Riko) (2) 

7V=0'^ 1=1 j=0 

where the operators R{p) diagonalize the problem on the infinite interval [— cxd, oo]. The function qn can be cast into 
different forms [s^l- One specific representation is given by 

..(M, {k}; X) = (-c)"exp[z.(Eo" fe.-Efp.)] (3) 

llm=l(P™ - km- tejiPm " Km-1 " «e) 

The above perturbative expansion is a quantum version of Rosales expansion (33j.(34|. Thus, 

J 2^^l4ije +c J 2^ (^,_^^_,,)(^3_^, + ie) + 

The inverse expression (direct Gelfand-Levitan transform) is also available (see [35| for the finite interval), 

i?t(^) = -^ / da;*^(x)e-*«^ + / dxidx2dx3g2{xuX2,X3;q)'f\xi)^\x2)'i'{x3) + ... (5) 



2tt J V2n 

where e.g. ^ 52(2:1, X2, xa; g) = c9{x2 - X3)9{x3 - xi) exp{-iq{xi + X2 - X3)), g4{xi,X2,X3,X4,X5]q) = c^9{x3 - 
X5)0{x5 — X2)0{x2 — X4)9(x4 — xi) exp(— i(7(a;i + X2 + X3 — X4 — 2:5)). Here the operators of refiection coefficient 
i?(f),i?^(f) diagonalize the Hamiltonian and satisfy the Zamolodchikov-Faddeev algebra (see Appendix A for more 
details) 

[H,R\0] = eR\0, (6) 

R{£,)R{0 = si^ - om')m), RmHa - ^(e - ORHORio + 2vr5(e - a- (7) 

where the scattering matrix is S{^ — ^') = (C — C ~ *c)/(^ — ^' + ic). These relations make the problem of finding time 
evolution easy: the factor exp[zx(^p ^j~Si Pi)] ^'i- © simply needs to be replaced by the factor exp[ix(^Q ki — 
'^1 Pi) — '''^iJ^o kf — Pi)]- Application of this formalism to non-equilibrium problems is also straightforward: 
one has to decompose the initial state, written in terms of bosonic \l/(a;)-operators, into a series of Zamolodchikov 
-Faddeev operators R{£,), generated by the inverse transform, and then find the time evolution according to ^ using 
the direct transform. 

Strictly speaking, this expansion is valid only for an infinite interval. For finite intervals it gives rise to singularities 
in expressions for the correlation functions, which, however, can be corrected by careful consideration of expressions 
in each order i3_2j. This series can be summed up explicitely only for the infinite value of the coupling c [37|.[3^. 
Unfortunately this approach is very difficult to use for calculating time dependence of the correlation functions at 
finite c even in equilibrium. We are not aware of the reproducibility of the asymptotic results which would correspond 
to the Luttinger liquid power-laws. On the other hand the advantage of this approach is the possibility to generalize 
it to other specific boundary or initial conditions. In principle this formalism should allow to include impurities 
and can be extended to multi-component generalizations of the NSE [39| . 
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B. Intertwining operator 

Some integrable (and sometime quasi-exactly solvable) models have the following property 

I Ho = HJ (8) 

where two different Hamiltonians i/o.c are connected (intertwined) by the action of some operator /, called the 
intertwining operator. In general, if Hq and He belong to two different representations of the same algebra, then 
intertwining operators exchange these two representations of the same algebras (or, saying mathematically, establish 
a certain homomorphism between them). If Hq is a Hamiltonian for free, noninteracting particles, and H^ is the 
Hamiltonian for the interacting system with interaction constant given by c, then the evolution of observables in 
the interacting model can be related to the evolution of the non-interacting one. Hence dynamics of the interacting 
Hamiltonian can be mapped to dynamics of the noninteracting Hamiltonian using 

e'"-^' = Ie'"°'r\ (9) 

Few comments are in order. The existence of the intertwining operator for integrable models follows from the 
following facts: (i) the Bethe states of all integrable models are parametrized by the integer numbers which label the 
eigenfunctions of noninteracting problem. Therefore the wavenumbers of interacting model are analytically connected 
to the wavenumbers of the noninteracting one; (ii) in the language of the coordinate Bethe ansatz a system of 
N particles is described by the wave function defined in the A'^-dimensional space divided by the hyperplanes on 
which the collision processes occur. Outside these hyperplanes a system behaves as noninteracting. Transition 
amplitudes between A^! different regions (outer space of the hyperplanes) are the same for all eigenstates because of 
the permutation symmetry. These transition amplitudes define a unitary transformation which is nothing but the 
intertwining operator. However, one should note that in general the operator / is not unitary. Thus, the wavefunctions 
of exactly solvable models are not orthonormalized, their overlaps are given by the determinants of some matrix via 
Gaudin-Slavnov formula Defining orthonormalized wavefunctions one can construct a unitary version of the 

operator /. 



1. Degenerate affine Hecke algebra: first quantized notations 

In the case of the Lieb-Liniger model the intertwining operator can be related to the representation of the degenerate 
affine Hecke algebra [i^]. For the XXZ spin chains there should exist similar intertwining operator between different 
representations of the Uq{sl2) and Temperley-Lieb algebras. One reason to expect this is because both LL model and 
XXX spin chain share essentially the same i?-matrix satisfying the Yang-Baxter equation. Some earlier discussion on 
the XXZ chain are contained in [4l| . 

The construction of / for the ID Bose gas relies on the notion of Dunkl operator 

d^ = -id, + i^ ^(e(a;, - Xj) - l)sij + i^ ^{e{xi - xj) + l)s,,j (10) 

j<i j>i 

where e(x) is a signature function, operator Sij provides a representation of the Artin's relations for the braids 
and exchanges coordinates of particles i and j and satisfy XiSij — SijXj. The operators di,Sj = Sjj+i form a 
representation of the degenerate affine Hecke algebra. Another representation of the same algebra is formed by the 
ordinary differential (difference) operator —idi and the integral operator Qi (4^ representing scattering matrix and 
acting on the arbitrary function /(. . . ,Xi, Xi+i, . . .) as 

Qif{. . .,Xi,x,+i, ...)=/(.. .,Xi+i,Xi, ...)- c f{...,Xi- t,Xi+i +t,.. .)dt (11) 

Jo 

The intertwining operator / interchanges these two representations of the affine Hecke algebra, [di, Sij) and {—idi, Qi) 
and thus intertwines the Dunkl operator di and the ordinary partial differential operator. Explicitly 

/= ^ 0{x^-Hl) < . . .X^-l(^N})Sw-^Qw (12) 
w&Sn 

where s^-i = s^^ . . . Sj^Si -, an d Qw — QiiQi2 ■ ■ ■ Qip (1 ^ *27 • ■ • , *p ^ N — 1) and where w is a transposition from 
the symmetric group Sm |40| . 



5 



All conserved quantities I„ for the ID Bose gas system are given by the powers of the Dunkl operator, In = ""('^D 
where 7r(.) is projection onto symmetric subspace. Thus, the Hamiltonian ([T]) is just equal to 

As a side remark we note that these facts are convenient for the formulation of the generalized Gibbs ensemble 
(GGE) approach to non-equilibrium dynamics of the LL model '55]. In particular, the GGE density matrix of the 
ID Bose gas should be related to the noninteracting density matrix by the intertwining relation. As discussed in 
Ref. I43i], the GGE conjecture should work when eigenvalues of the integrals of motion are parametrized by either 
{0, 1} (fermionic-like systems) or by integers {0, 1,2,.. .} (bosonic-like systems). The case of the Bose gas belongs to 
the second class. We therefore conjecture that the GGE is applicable to the ID Bose gas (LL model), at least for the 
local observables. 



2. Graphical representation: Gutkin approach 

An interesting approach to the ID Bose gas based on the intertwining operator has been developed by E. Gutkin in 
a series of papers [4l| , [1^ , [3| , [i^ . It was later applied to the problem of propagation of an optical pulse in nonlinear 
media in Ref. (46| . 

In this approach the intertwining operator is expanded as a series of different contributions labeled by the so-called 
collision graphs ^9^. The evolution of the quantum field ^(a;,t) governed by the NSE is given entirely in terms of 
the evolution of the field 'i/o{x,t) for the free Schrodingcr equation, idt'^Q{x,t) = —d^'i>o{x,t), and by the quantities 
Or (defined in terms of some distributions) entering the expansion of the intertwining operator over a set of collision 
graphs 

M"a; M"yar(xi,X2,--- ,x„,yi,y2r--yrO*S(a;i)*S(x2)---*S(x„)*o(j/i)*o(j/2)---^'o(y«) (13) 
r,n=g(r)"^ 

and similarly for the operator with ap replaced by the related distribution br- There is a special choice of or and 
br for which the operator I is unitary. For details about this approach we refer to the original papers by Gutkin. 

The advantage of this approach is the possibility to find the propagator of the nonlinear problem and thus, poten- 
tially, solve the initial state evolution problem for an arbitrary initial state. Computationally one only needs to deal 
with the free fields i^) and their time evolution with a free (noninteracting) Hamiltonian. We note that the collision 
graph expansion is highly non-perturbative and can not be rewritten as an expansion in powers of c in opposition to 
the Gelfand-Levitan-based expansion of Sec. Ill Al A disadvantage of this approach is rapidly growing complexity of 
collision graph expansion with increasing number of particles. Also it is not known at present whether some of the 
most relevant graphs can be summed up in order to develop useful approximation schemes. 



3. Second quantized approach 

In series of papers [13] Sasaki and Kebukawa developed a field-theoretical approach to the LL model. Although 
they did not present it in this context, the second quantized form of the intertwining operator can be recognized in 
their construction [97|. In this approach one considers the LL Hamiltonian written in the usual second-quantized 
form, 

2 ^ 
P P:Qjr 

Here ajj,ap are bosonic creation and annihilation operators corresponding to momenta pi = 2TT'hni/L [ni should be 
integer). The second-quantized analog of the Bethe ansatz wavefunction for N particles has the form 

TV 

i*?i,92, ...?«)= 59i,g2,...9iv X! n ^(^^j'^^j^n^i:" 1 ■ ^+9='^'' 

Pi.j;l<i<j<N l<i<j<N i=l '-^-'^^ 

where Bq-^^q^^ ^ q^^ is a normalization factor, pij = Pi—pj = —Pj,i and (i = 1, . . . , N) are defined as an integerx27r?;,/L, 
d(j)i^j] ki^j) = —kij/{pi_j — kij) and kij = —kj^i are solutions of the Bethe ansatz equations 

cot{^^)^—{k,-k,) = —{2k,^j + y2{h^i-kj^i) + q,-q,), l<i<j<N. (16) 
2Ti VIC mc ^ 
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The eigenvalues in these notations are given by = J2i=i kf/2m = J2i=i J2j^i{ki,j + 9i)^/2m. 

It is interesting that authors of [J?] found a unitary operator which transforms eigenstates of noninteracting system 
Hq into eigenstates of the interacting Hamiltonian H^. In the second quantized form this unitary operator is given by 

n n 

qi^q2^---qn Pi.i\^^i<j^'n l<i<j<n i—\ i—1 

where are normahzation factors for non-interacting eigenstates. ExpUcit calculations show that this unitary 

operator is indeed equivalent to the intertwining operator from previous subsections. 

This approach has advantages for initial states which can be readily represented using formalism of second quanti- 
zation. For the case of a large coupling constant it gives the same result as the one used later on in the text. 



C. q-bosons 

In a series of recent papers [11], [1^ a system of g-bosons hopping on ID lattice has been studied. This model is 
defined by its Hamiltonian 

1 ^ 

Hq^--Y.{hih,,+i+hX+i-'^N^) (18) 

where the periodic boundary condition is imposed. Operators _B„, B\ and iV„ satisfy the g-boson algebra 

[N,,b]]^b]6,j, [N,,b,]^-bA,, [h,b\] ^q-^^^S,, (19) 

where q = . The operator of the total number of particles N = X^^li commutes with the Hamiltonian Hq. 
The continuum limit of the model is defined using the limiting procedure 

c5 

(5^0, MS^L, 7=Y- (20) 
When applied to the system of g-bosons this procedure gives 

H=[ dx [d^b'' {x)d^b{x) + cb'' {xp {x)b{x)b{x)] , (21) 
Jo 

where [b{x),b^ (y)] — 6(x — y) are the canonical Bose fields. Therefore ID Bose gas with contact interaction can be 
regarded as a continuum limit of the g-boson lattice model. Apparently a solution of the inverse scattering problem 
for the NSE should be related to the one for the q-boson model. An explicit form of this relationship is not known, 
although should exist [osj . 

This approach can be used to study non-equilibrium dynamics. Solutions of particular evolution problem for q- 
bosons can be transformed into solutions of NS problem via the limiting procedure l|20p . An attempt to construct 
the evolution operator for a single-mode g-boson system has been made in [s^]. Even for this simple case the g-path 
integral has an involved structure which includes integration with the measure corresponding to non-flat phase space 
and complicated action having non-trivial dynamical phase. All this makes the possibility of explicit evaluation of 
the g-path integral questionable. It is interesting to note that dynamics in this phase space is naturally constrained 
by integrals of motion. This is reminiscent of recent suggestions of the importance of GGE for dynamics of integrable 
models 221. 



D. Extended conformal symmetry 



It is known that the low-energy description of the ID Bose gas can be formulated as a bosonic Luttinger liquid 
[Sll [52j . This description operates with a linear spectrum of bosonic modes and treats the interaction part essentially 
exactly. From the more formal point of view, in this low-energy description the thermodynamic limit of the NS 
system represents a conformal field theory (GET) with the central charge c = 1. This means, in particular, that the 
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low-energy (bosonic Luttinger) Hamiltonian is a linear combination of the zero- momentum generators Lq^Lq (for the 
right and for the left moving parts) of the conformal symmetry algebra, i.e. the Virasoro algebra 

[L„, Lm] = {n- m)L.n+m + ]^"-("-^ ^ l)5n+m,,Q (22) 

and the same for L„, ([-L„, Lm] — 0, Vn, m). The space of states is the so-called Verma module (representation space) 
for this Virasoro algebra. The Kac table defines a conformal tower for this c—\ CFT which allows to compute any 
correlation functions. One could also include finite size corrections and effects of irrelevant perturbations [ssj.fs^. 

The non-equilibrium dynamics can not be correctly described in terms of the low-energy, Luttinger liquid theory 
only. The reason is two-fold: first, if the initial state represents a "large" perturbation over equilibrium, which includes 
the overlap with the whole spectrum, the excitations involving nonlinear part of the spectrum should be important. 
Second, even for relatively weak perturbations, which do not thermalize because of integrability, the dynamics at large 
times leads to growing correlations over the entire momentum space. So irrelevant contributions from the interaction 
part of the Hamiltonian should become important at long times. 

It appears that nonlinearity of the dispersion as well as irrelevant (from the equilibrium point of view) parts of 
interaction can be included using the extended conformal symmetry called Wi+oo- The Wi+oo algebra is a represen- 
tative of a class of nonlinear algebras which appeared in conformal field theories. It is a product of a quantum version 
of an algebra of area-preserving diffeomorphisms of a 2D cylinder and the abelian Kac-Moody algebra. The presence 
of dynamical symmetry Wi+oo is a general property of gapless one-dimensional systems. The chiral generators 
of Wi+oo are labeled by the momentum index n, k — 2-11x1/ L and by the conformal spin index a = 0, 1, 2, ... 00. They 
satisfy 

[W„", W^] = {Pn - am)W:+t^ + g(a, /3, n, m) W,^+fr' + • • • + 5'^^ 5.n+rn^^cA(a, n) (23) 

where q(a, /3, m) and d(a, n) are known polynomials of their arguments [ssj . 

For many interesting cases the Hamiltonian can be written as a linear and bilinear combinations of Cartan generators 
of the Wi+oo algebra 99J. The representation theory of in the Verma module can be constructed in the same 

way as representation theory for the Virasoro algebra [s^ starting from the highest weight state and building a towers 
of descendant states. Using this representation theory one can construct the action of Wi+00 generators on bosonic 
Fock space. One can therefore establish a correspondence between Bethe ansatz eigenstates and certain combination 
of states of representation module of the Wi+oo- The evolution operator is factorized in the product of factors 
corresponding to different Cartan generators 

U{t) = e'"' = \[e'^-'^''' (24) 

where Pa are functions of the interaction strength and a level. When the initial state can be expressed as a linear 
combination of descendant states, the application of U (t) is straightforward. 

In practice it is difficult to deal with the whole nonlinear dispersion, so, for practical purposes one can keep only 
next-to- linear terms in this expansion. In this case the Hamiltonian will include only a finite number of Woo generators. 
For example, in the strong coupling regime, corrections to the Hamiltonian to the first order in the inverse power of 
the interaction strength, 

^ = + ^« ) + +1^0^] + ... (25) 

The advantage of this approach is simplicity of including higher-order corrections originating from nonlinearity of 
dispersion, and the potential to solve the problem exactly. The disadvantage of this method comes from its limitation 
to treat gapless systems only. Details of this approach together with several related questions will be presented in a 
separate work |57| . 



E. Form factors and decomposition of the initial state 



This is the most direct approach which we will use further in this paper. In this approach we decompose the initial 
state in terms of the eigenstates of the Hamiltonian. The latter form a complete and orthogonal set of states. To 
compute time evolution of the correlation functions we have to combine several ingredients: (i) find complete basis of 
many body wave functions; (ii) know exact eigenvalues; (iii) determine matrix elements (or form-factors) of various 
operators in the basis of exact eigenfunctions; (iv) find decomposition of the initial state in the complete set of exact 
many-body states; (v) develop effective procedure for summation over intermediate states. 
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Ingredients (i)-(ii) are known from the Bethe ansatz exact solution. Matrix elements (iii) were computed in many 
exactly soluble models on the basis of the so-called determinant representation. Recent progress in computation of 
the matrix elements of operators in the basis of the Bethe Ansatz wave functions makes it possible to advance in 
this direction. Decomposition of the initial states (iv) over the complete basis depends significantly on the concrete 
nature of the states. We were able to evaluate these overlaps for a ID Bose gas at large interaction strength for 
two types of initial conditions: all particles in a zero momentum state and all particles in Gaussian wavefunction in 
real space. Currently we perform summation over intermediate states numerically. This imposes certain restrictions 
on the number of particles in our system. We checked that in equilibrium this number is sufficient to saturate any 
correlators to their values in the thermodynamic limit. To illustrate this observation we plot the g2{x), computed in 
the ground state using our method, and compare it to the known analytic expression in Fig.([T]). 

Here we make general comments on the structure of the phase space of form-factors. For a systems with a gap 
only a small number of states needs to be taken into account. Contributions of many-particle states are suppressed 
by the gap. This is the case for the NLS system with attraction, which can be identified with the nonrelativistic 
limit of the quantum sine-Gordon model for which previous statement is also correct (see e.g. review [i^l and refs. 
therein for direct comparison of different contribution to the form- factor expansion of the correlation functions). 
Contributions coming from higher order terms grow in the limit of long time evolution, but should only generate 
subleading corrections. The situation is very different for gapless systems, such as the NLS system with repulsion. 
Many-particle contributions are suppressed at most as a power law. Hence ideally we need to take into account the 
entire phase space. This makes the problem of summation over intermediate states very difficult. In this paper we deal 
with the most difficult case of repulsively interacting Bose gas in the regime of strong interaction, when, in principle, 
the allowed phase space is huge and multi-particle states are in general not suppressed. 

The form-factor approach has been successfully applied to computation of equilibrium correlation functions in many 
models [slj. A review of the applications of this method to massive models is given in Ref. [H^. The gap in the 
spectrum leads to the rapid convergence of the form-factor expansion. It is enough to take contributions from a few 
particles to saturate any correlation function. This is not the case for massless models where contributions of multi- 
particle processes are essential. Recently this approach was applied to compute equilibrium correlation functions of 
ID Bose gas in [13], [HI, (see also Ref. ^] for the attractive case 100]). Earlier studies reviewed in jjj allowed to 
evaluate various asymptotics of the time-dependent correlation functions at equilibrium (for concrete example of non- 
zero temperature see e.g. (63j. l63 .) Recently, a form-factor approach has been applied to time evolution of specific 
initial state in the XXZ model [65| . Specific choice of the initial state made it possible to perform analytical analysis 
of the time-evolved. This is one of the rare examples when a work function and Loschmidt echo, important quantities 
related to the properties of a stationary state [6a], have been evaluated analytically. 

One of the advantage of the direct application of the form-factor technique is the possibility to consider non- 
equilibrium dynamics of weakly non-integrable systems. In Refs. [67| the form-factor perturbation theory has been 
developed for a class of models which deviate weakly from an integrable " fixed point" . This perturbation theory is 
unusual because its unperturbed states are highly correlated states of interacting integrable theory. We expect that 
this formalism can be extended to treat non-equilibrium dynamics as well. 

The form-factor approach supplemented by a procedure of decomposing the initial state into a complete set of 
many-body eigenstates is universal and can be applied to a large class of integrable models. This method allows 
several extensions and modifications. In the rest of this paper we focus on this method and demonstrate that it is a 
powerful tool for analyzing non-equilibrium evolution of correlation functions. 

In the next section we also show that one can use the intertwining operator to find an ov erlap coefficient between 
interacting and noninteracting states. We do this explicitly in the strong coupling limit. |lOl| To conclude this section 
we note that other approaches to calculating correlation functions in integrable systems undergoing nonequilibrium 
dynamics have been discussed in Refs. [6^, and [g^. When initial states can be written using ordering of the 
coordinates as in eq. (|106|) (see Appendix A), an approach to calculating their dynamics has been proposed in Ref. 
[to} . We expect that the approach developed in [t^ is related to the method which we discuss in the current paper. Yet 
another approach, closely related to the form-factor method in the limit of an infinite volume, allows for computation 
of the time evolution of the one-point function. This method relies on representing certain (integrable) initial states 
as boundary states, which correspond to coherent superpositions of the eigenstates of the model [7l|, [72l. Recently, an 
interesting approach to equilibrium correlators in the Lieb-Liniger model has been developed in Ref. [73j ] . The idea of 
the latter papers is to consider first a relativistic massive model (the so-called sinh-Gordon model), which, in specific 
non-relativistic limit, is equivalent to the Lieb-Liniger model. The form-factors of the sinh-Gordon model are known 
and the form- factor expansion of correlation functions is rapidly convergent. Such convergence is a result of the mass 
gap in the spectrum. In the special non-relativistic limit, all expressions for the field operators and correlators of 
the sinh-Gordon model go over into corresponding expressions for the LL model. Moreover, the non-relativistic limit 
benefits from the rapid convergence of the form-factor expansions in its relativistic counterpart. One might hope that 
this interesting approach can be extended to non-equilibrium dynamics. 
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III. EXACT APPROACH TO NON-EQUILIBRIUM DYNAMICS BASED ON FORM-FACTORS 

In this section we describe the apphcation of the method of form-factors to the computation of non-equihbrium 
correlation functions. This method is general for any integrable system. In this paper we focus on the ID Bose gas. As 
discussed earlier the treatment of gapless system is more challenging because of the large number of contributions from 
intermediate states which need to be included. To overcome this difficulty we devised a special numerical procedure 
of summation over intermediate states. This procedure is described in the next section. 



A. Formulation of the problem 

To fix notations we consider the nonlinear Schrodinger Hamiltonian on a finite interval [0,L] 

H=[ + c^'t(a:)^'t(j;)*(j;)*(j;)] (26) 

Jo 

where the commutation relation between bosonic field ^'(x) and its conjugated 4'^(x) as well as the pseudovacuum 
state |0) are defined as 4] 

[^{x),¥{y)] = S{x-y), [^{x),^{y)]^[^\x),^^{y)]^0 (27) 
^'(x)|0) = 0, (0|*^ = 0, (0|0) = 1. (28) 

The total number of particles and the momentum operators are conserved quantities 

[ *1"(x)*(2;), P^-l-j [^\x)d.^-^ -{d-^\x))-^{x)] (29) 
[H,N] = [H,P]=0 (30) 

In each iV-particle sector this system is equivalent to the ID Bose gas with contact interaction potential (Lieb- 
Liniger model) [l|, In connection to the previous section we note that inverse scattering transform for this 

problem has been developed in many papers (see [sl], 74[,[33l, 44 [ for extensive reviews). 

The problem we are asking here can be formulated as follows: suppose we prepared a system in a certain initial 
A''-particle state IV'o^'') (e.g. coherent state can be considered as a superposition of states with different N). This 
state evolves according to the Hamiltonian ((26|) . The question is to compute the correlation functions of the following 
type 

{4''^\'f^{x,,h)^>{x2,h)\4''^), (31) 
{4>^J'^\^Hxuh)^ix,,h)^Hx2,h)^ix2,h)\4''^), (32) 
(Vr^|exp(aQ(x))|4^)), (33) 
where the last formula expresses a generation function for the density-density correlator, 



Q{x,t)= / ^\y,t)^{y,t)dy (34) 
Jo 

The density is defined as 

p{x,t) = -i>^{x,t)^{x,t) (35) 
Assuming that the initial state is normalized, one can insert a resolution of unity several times, 

into expression for the correlation functions. This suggests to define normalized correlation functions, 

{4''^\0{x,,t,)0(x^,t^Mr) 

= V V V (^r^l{^}A-)({^}^|0(^i'^i)liA^}A-)(M^|Q(^2'^2)|{^}Ar)({4jv|4'^^) , . 

m i^iTT ({a}|{a})(MIM)(MIM) ■ ^ ' 

tA}]v iM/iv {i^Ijv 
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In the following we concentrate on the density-density correlation function. According to our programme we need 
matrix elements of the operators in the complete basis of states. In this basis time and space dependence are trivial: 
using the Galilei invariance (we implicitly assume periodic boundary conditions) the matrix elements which we need 
are given by the expressions 

mN\Q{x,t)\{X}^} = (e-(^r'-Pr') - l)e^'(^r-4"VQ({M};{A}). (38) 
Here the energy and momentum are given by 

N N 

i^f^^E^^ pr-T.^. (39) 

In the following it will be convenient to use more compact notations 

4°^ = (V^r^KA}^) (40) 

^q(M;{4) - (M^IQ(o,o)|{4jv) (41) 

ll{A}|| = ({A}|{A}). (42) 
We can write down an expression for the density-density correlator as 

(4~V(^,i)p(0,0)|^r) (43) 
^ d' V V V 4°^(-4i°VfQ({A};M)(Fg(M;{^}))* 
d^^dx.'f- f-^ ll{A}||||MIIIIMII 

|A>jv IM/N \V\n 

X (e-i(Pr'-^r') _ i)e^ti(4"'-<')(e-^(^i"'-^r') - l)e^*^(^^<"'-^=r') (44) 
where * means complex conjugation. 



B. Bethe Ansatz Ingredients for ID Bose gas 

As discussed above in Section pi El) we need several ingredients: a set of exact eigenstates and eigenenergies of the 
model; an overlap of the initial state with the eigenstates; an explicit expressions for the matrix elements of operators 
in the eigenbasis. 

1. BA states 

The BA states are described by a set of A'^ real numbers A which are given by the solution of the BA equations 

I A, - A, 27r , ^ A + 1 , 



A, + -}_^2arctan^ L ^ _(/^. _ __), j^l...N (45) 

1=1 ^ 

Here the quantum numbers Ij are half-odd integers for N even, and integers for N odd. For eigenfunctions rapidities 
do not coincide, Xj ^ Afc for j ^ k. The whole Fock space is obtained by choosing sets of ordered quantum numbers 
Ij > Ik, j > k meaning that Xj > Xk, j > k. 

From solution of these equations we immediately obtain energies and momenta via Eq. (j39|) . 

2. Overlaps between BA states 
The overlap between BA states is given by the Gaudin-Korepin formula (zll, (z^ Q, [i,[z3i [zi|,[zi|,[ili,[6li 

({AV|{A}^.)=c^ n %r^^r^^^^^«^» (46) 
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where 

N 



1=1 



K{X„Xk) (47) 



^ (A— ^''^ 



3. Matrix elements 

The matrix elements (form-factors) for the current operator are given by the following determinantal expression 

MM. M 



^q(m.;{a}.) = n ( '^:\^" ) 



X det SM- -^) + ,r ~ ^'^^ - ^(^^ - ^'^^J (49) 

where 

^ nf^i(A^a-A, ±»c) ^^P^ 

na=l(Aa - Xj±tc) 

and where if is given by p8)) . The matrix inside det() has only real entries. This matrix has size N x N. Real 
numbers {A}, are solutions of BA equations pS]). 

C. Overlaps of the initial states with Bethe Ansatz basis 

Here we study in details two physically motivated examples of initial states. The first one has its origin in the 
condensate physics whereas the second one describes a Gaussian pulse created in special nonlinear media. Before 
going to these examples we make several general observations. 

1. General remarks 

To compute the overlap with the initial state we consider the following (coordinate) representation for the Bethe 
ansatz states of the NSE 

|*Ar(Ai, . . . Ajv)) = / d^zxNizi, . . . , zn\Xi: . . . , Aw)^'t(zi) . . . *t(^^)|o) (51) 



where the function xn is given by 

d d 

Xn{zi,. .. ,zn\Xi, ■■■ ,Xn) ^ const Y[ (^ - ^ + c) det[exp(iAjZfe)] 

N 



: ^(-l)'^l exp[i ^ ZnXvn] W[Xvj ^ Apfe - ice{zj - Zk)]. 



n,>;c [(A, -xuY+ V :ri 

(52) 

It is convenient to define the basic wave function 

|vI/)(^) = _L.vI;t(^^)...vI;t(^^)|0). (53) 
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In principle, any quantum state can be expanded as 

1*^°^) =E«" / -^fn{xu.-.,xMi"^ (54) 
„ J vn! 

with 

^|a„p = l, / |/„(a;i,...,x„)pdxi ...dx„ = 1 (55) 

Intuitively one expects that the evolution of some special initial states can, under certain assumptions about these 
states and about the Hamiltonian, be relatively simple. This is indeed the case for the ID Bose gas. Using the 
expansion ([2]) and contour integration it was shown in fs?] that the state defined as 

\i^N{t = Q))ord = 0{xi > X2 > . . . > xn)\^\xi)¥{x2) ■ ■ ■ *t(a;w)|0) (56) 

is equal to the state R'' {xi)R'' {X2) ■ ■ ■ R'' {xn)\0) where R{x) = J ||e"«i?(C) (see Section HEl for definitions). The 
calculation of the evolution of this state is therefore straightforward. 



2. The overlap with a delta function in momentum space. 



|*(gi),...vl/(g^)) = — i=[/ e'«i"iM't(xi)dxi],...[/ e'''-^-^\xN)dxNm. (57) 



We first construct a Fourier transform of the basic state: 
At the end of all computations one could take all momenta to be equal qi = q2 = . . . qn io obtain a condensate state 

\4''\'i))c = [^Hq)f\o)- (58) 

Let us consider the overlap of this state with the Tonks-Girardeau wave function corresponding to the large c 
limit. There are several reasons to focus on this limit: 1) from our basic interest in application to quantum optics 
with strongly-correlated photons the TG case is the most interesting since effects of fermionizations ^] should be 
most pronounced 20] in this limit; 2) equations for the overlap are relatively simple and amenable to analytical 
calculations; 3) one could use an expansion in powers of 1/c around TG limit to produce results for the correlation 
functions beyond TG limit. 

In the TG limit the wave function takes the following form jl^, [83j.[83|: 

1*^) = i r dx, ... r dxNX^TGi{^^}\{^^})'^\^^) ' ' ' "^^^ nW) (59) 

vJyi Jo Jo 

such that ~ 1 and where 

X^G^({MIM) = ^=L= n <^.-^.-)det(e"'"^") (60) 

Here 

T"" \ e Z + i, N even 
The overlap therefore is reduced to the following integral: 

cL 



2tt Jn, eZ, N odd; 



c{4''\q)\^TG)^j^ I dxi I dx2... I dxN n e(2:,-Xfe)det(e'^"(^"-«)) (62) 

■^0 Jo Jo l<j<k<N ^ 

where the factor iV! comes from the overlap with the basic state and denominator comes from the normalization. 
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The result for the overlap is given by 

c{4''\q)\^TG) = (63) 
^ Nl2^Ui<,<k<Ni>^j-^k) ^ j (_j)f(cos(^))f(sin(^))f, N even; 

i^nL(A.-g)ni<.<,<^(A. + A, -29) 1 (-z)^(cos(f ))^(-sin(^))^, N odd 

X (_l)i:f..n,g^p(!_^^^.)g^p(_^^^|) (g4) 
J = l 

This formula can be simplified further — see Sec. IIV C2\ below. 

Straightforward generalization is a construction of a coherent state in fc-space. One can take a superposition of 
I'^'^^Hl)) states with different iV's for real or complex a: 

00 TV 

l«(9))-e-^l"l^^^|*W(g)) (65) 

We will not consider this state here. 

3. The overlap with the Gaussian pulse 



In the case of a pulse prepared at time t = we consider a state 

Xl,. . .XN 

)*t(^i)...^,t (2.^)10) (66) 

i=i 

and choose a function f{xi...XN) in the form of a product of gaussians, 

N 

/(a;i,...xjv) = nexp[-a:?/2a2] (67) 

1=1 

For simplicity, we take cr,; = ct, V«. The overlap with the TG state is given by 

(*TdV.r>G = ^(i)"/^-"exp(-^5:A?)n[Erf(^^^^) + zErfi(^)] (68) 
For all practical purposes the Erf function can be replaced by the gaussian. 

D. Tonks-Girardeau limit 

From our understanding of the equilibrium structure factor of the ID Bose gas [i^ we expect that the allowed phase 
space of the weakly interacting gas is limited to the vicinity of the dispersion curve of Bogoliubov quasiparticles. On 
the other hand we expect that in the opposite case of the TG limit the allowed phase space should be considerably 
larger. Hence non-equilibrium dynamics should be more involved in this case as well. Here we thus focus on this 
interesting limit in more details. Certain simplifications occur in the equations for the matrix elements, which are 
described below. 

1. Reduction of the form factors 
To compute density-density function in the TG limit 

{4^^\pixi,h)p{x2,h)\4'^^) = -^5^ J2 FP;{xiM,\lJ)FP,{x2M,l^,y)AxA, (69) 
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we have to know the following form-factor: 

i^j^(x,i,{A},M) = e'*(^^-^-)— ^ / dxi...dx^^(-l)[^l+[21cxpH^(Ap-A^Q)] (70) 



N-1 



exp[-i ^ Xa{Xva - AisJ]- (71) 







A limiting procedure for the form factors in this case has been developed earlier [85[ , [86| , [87[ . We use here the fact 
that the TG limit can be obtained as a double scaling limit of the XX chain. Using 

'''^dxe-"(^»-''^) = I f" ^ = ^^ (72) 

[ 0, Aa Mb 

where we use the boundary conditions e'-^'**" — e'-^^*" = (—1)^+^. Therefore 

F^ix,t,{X},{^^}) = \ f/X-^.ix-,) fl^,/''^ X n ^L^{p,5x,+ '—^{l-5^,)) (73) 

\^ L e ^ Ai — ^1, . . . Aw-i — ^jv-i, A 7t ^ I, 

and otherwise. 

Note that this form-factor is antisymmetric function with respect to any interchange of momenta for both sets {A} 
and {fi}. So when they are not ordered in the formula above, we have to assign a factor {-l)[P]^+[PU ^ where [P] IS a 
permutation index. 

2. Delta function subtraction 

When we calculate density-density correlator (|69| . we need to subtract a contribution of the particle with itself. 
This contribution comes from commutation relations for Bose operators: 

{X\¥{xi)-^\x2)'i'{x2)'f{x,)\l^) - {X\¥{xi)'if{xi)^\x2)'f{x2)W) - {X\'^Hxi)'^ 5{xi - X2) (74) 

Considering the TG limit and using (|73p we can calculate the corresponding (5-function contribution 

- {X\^Hxi)'${x,)\j,)S{xi~X2). (75) 
For the ground state in the continuous limit this contribution is given by 

-^J dqe^^'-^'^^ (76) 

where x = Xi — X2 and t = t2 — ti. In general case (169^ delta function yields an additional term 



\po + \ (-l)[^l^+[^l-AAA,e'-^(^"-''^)-^*^(^--^"M ^5]e"'- 



itE„ 



(77) 



Now we can combine all components together to calculate the density-density correlation function in the TG limit. 
For the ground state we have 

g^^\Ax,At)^p^ + -^ f f e^*(??-«2^)cos[(<zi-<Z2)^]-;l r"dge*«'*-^«^ (78) 



^" ■'\qi\>-lTp J\q2\<Trp ^" J -oo 



In a general nonequilibrium case correlation function can be written as 

{^Plf\qMx,,t,)pix2,t2M^^\q)) =p'o + ^ E |A,pe^-(^"--^)+''*(^-^^) 



^P^— ^ (_2)[-P]A + [P]M^^y^^(g^^l(^a-Aib)-lil{£^A-£^M) ^1X2 {Xa - fJ^b) -'tt2 (E x - E ^ 

a.b-.Xay^f-i'by^ 

-[po + 2 E (-l)[^l^+[^l-^AA^e*"^(^«-^^)-**^(^^-^-M -^E^'^ 



^ixq — itEq 



(79) 
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where the first term represents S — S contribution, the second one corresponds to X = v (diagonal part), third and 
fourth terms come from the nondiagonal {X ^ fj, ^ v) part, and the last term is a delta contribution discussed above. 
Here everywhere A differs from fj, by one filling number only, and fi differs from v by one filling number as well. 

Here, the S-S part as well as diagonal parts are symmetric with respect to interchange of momenta since these are 
the products of two antisymmetric functions. On the other hand in off-diagonal terms we inserted the sign factors 
in order to ensure that proper symmetry is preserved: since each F is antisymmetric, the nondiagonal part will be 
antisymmetric as well provided that functions A\ and A^, are symmetric. 

The second term (diagonal part) can be simplified further using the following relation 

N N 

E /(A.)=5^/(A,)-^/(m,) (80) 

-^aT^Ml ■■■Mjv a—1 b—1 

for arbitrary function /. Then 

1 1 ^ 

|v4Ape*^(^"-''^)+^*(^--^-) =--^|^AAe'"^''-"^-|2 (81) 

a,b:Xa^l^i,yb a=l 

-. N N 

a=l 6=1 

This identity can be also used to simplify the nondiagonal part. 



3. Deviation from the TG limit 



The main difficulty of our approach is computing overlaps with the initial state in a compact, analytical form 
for arbitrary interaction strength. At present we do not have a complete solution of this problem for an arbitrary 
state. For particular types of states the problem can be attacked using one of the formalisms of Section II. Thus, 
in particular, if the initial state can be expressed in second-quantized notations via creation-annihilation operators 
acting in momentum space, then the second-quantized version of the intertwining operator can be used. If the initial 
state has a simple form in coordinate space, the coordinate version of the intertwining operator might be useful. One 
form of this operator in real space, mostly useful for the 1/c expansion around TG limit, has the form of differential 
operator acting on the TG wavefunction 0,(1^. It relates the TG wavefunction (which is almost fermionic) with the 
finite-c Bethe state 

^ l<i<]<N ■> 

Manipulations using intertwining operators will be discussed in more details in the future in connection to c oncrete 
physical problems. Here we note that expanding the intertwining operator up to order 1/c one observes |l02j that a 
correction to the overlap coefficients up to the order 1/c (neglecting the boundary effects) for the condensate initial 
state have the following form 

4°) ^ Af\l + !i^^pf)) + o(i/c2) (84) 

Note also that deviations from the TG limit includes other ingredients. One of them is an expansion of rapidities in 
powers of 1/c, 

A, = (2j-7V-l)^ 

Using this expression one can get a systematic expansion of the form-factors. Wc do not need these 1/c corrections 
since our procedure uses numerical expressions of the form-factors. 

We finally note that a computation of the overlap of the arbitrary BA state (0|C(Ai) . . . C(A„) with the initial state 
of a special form ^'(a;i) . . . 5'(a;„)|0) can be considered as a particular case of evaluation of the multiple point form 
factor if we regard the pseudo vacuum state as the BA state. The problem of computation of such a multi point 
formfactors can be solved using the formalism of multi-site generalization of the NLS problem. This problem was 
addressed in [69j . Another way to go beyond the TG limit is to use different formalisms described in Sec. II. 



4(^ 



' c 



A-K^ 1 

f + —N{2p + {N + l){N- mi^f 



(85) 
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IV. NUMERICAL TREATMENT 
A. Introduction 

In this part we present results of numerical computations of correlation functions for various non-equilibrium initial 
states and their time evolution. In order to do that, we needed to sum expression (1441) over all initial A, intermediate 
II, and final v states. 

Although we are interested in thermodynamic limit L — ^ cx), iV — >■ oo, numerically one can compute correlation 
functions for finite size systems only. In order to infer thermodynamic limit from our finite size results we can set 
the size of the problem {N and L) to be sufficiently large so that numerically computed correlation functions would 
be sufficiently close to the ones for an infinite system. One of the criteria for such sufficiency could be comparison 
of the numeric results to analytic results for the system for which thermodynamic solution for correlation functions 
is exactly known. For example that we already know the exact expression for a correlation function in the ground 
state ([75)1 . So we can use it to check the accuracy of our numeric procedures. 

Let's estimate the cost of a brute force approach towards summation (|44l) for a reasonably large system size. Because 
we are interested in non-equilibrium states, states in A, fi, and v are filled way above the "Fermi momentum" . We can 
estimate the number of terms in the sum (j44p by taking N ~ 100 and momentum cutoff ~ lOA^. Then the number 

of terms in the sum (^^^) - 10^20. Obviously, the sum can't be taken in a straightforward fashion. 

Now we would like to look for possible simplifications of the problem. The first hint comes from , where authors 
claim that when one considers summation of the form 

J2\{GS\F\i,)f, (86) 

where GS is a ground state, then the sum has a very limited number of major contributions, most of which are one 
particle. Note that this statement is not directly applicable in our case because we are not working with a ground 
state and also because we have a double summation, which might have different major contributions due to larger 
phase space, or not to have such major contributions at all. 

Second simplification might come from the fact that we are interested in a large interaction constant c case. For our 
purposes we can consider the case of large c. As we found in the previous section, relevant quantities become much 
simpler. The most prominent feature of the TG limit is that the only non-zero matrix elements (AI-Fj/x) come from 
zero and one particle processes (j73|) . It means that in this limit for a given A we need to do a two particle summation 
in (j69p with one particle transition A — > /i and one particle transition fj, ^ v. This observation also tells us that for 
large, but finite c, processes with small number of transitions should dominate. 

Finally, structure of overlap factors A\ for different initial states might provide some clues on how to simplify our 
task. In particular for the initial states we are going to investigate, overlap factors can be simplified analytically so 
that we do not need to perform the summation (j52p over all permutations. 

B. Implementation 

Now we would like to discuss ideas behind concrete implementation of our calculations. Almost everywhere we 
consider the limit of large c, though we can do numeric calculations for arbitrary c at a cost of substantially higher 
cpu time. 

For a ground state we do several large-c calculations to check what values of c can be treated as "infinite" and to 
investigate multiparticle contributions. For large, but finite c, it is possible to further simplify expressions (|46p - (|50|) 
by doing 1/c expansion, as explained in the previous Section. We implemented both direct and simplified versions of 
matrix elements and overlap calculations and did not find any visible discrepancy in the results, though performance 
improvement was substantial. 

After we sum up the expression 

^ AMMF\t^){f^\FW) (87) 



over intermediate states /i and final states for a given initial state A, we are left with a problem of summation over A. 
But the phase space is huge, therefore we can't use a direct summation. We use Monte-Carlo summation instead [s^. 
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Namely, we sample set of states {A'} from the set of all possible initial states {A} with probability of state A' to be 
selected being P\. Then we can replace summation in (j87p over A with summation over A' to get 



Ex'.,,.AyA,{X'\FUl){fl\FW)/Py 



As we increase the size of sample {A'}, the value of this sum converges to the true value of (j87l) for any probability 
distribution P\. Depending on the pick of the probability distribution function convergence time can be very different. 
For instance if we take a uniform distribution over the entire set {A}, than we can see from Eq. (j68p that weights A\ 
can be very different. As a result, many contributions to the (1881) will be negligibly small. On the other hand if we 
are picking A such that \A\\ is around its maximum, all the contribution to the sum will be substantial. 

Indeed, it is known [83, that the sum optimally converges to the true value if the distribution function is proportional 
to the expression we sum. In our case it would (roughly speaking) mean P\ oc A\ A^, where v differs from A at 
most by two particle process. As a reasonable approximation, we use Pa A\. 

The next question is how to sample such states A'. We cannot sample states from the distribution directly, but we 
can use Gibbs sampling (90| to do that approximately. Note that any overlap factor can be expressed as a function of 
state A and parameters of the problem. A itself is a function of a vector of quantum numbers / (|^5t . We can fix all the 
components of / except one (/j ). By varying this component and resolving the equation (P5|) relative to A, we can find 
conditional probabilities P{Ij\Ii . . . Ij-ilj+i . . . In) oc A^^^j , where X\Ij is a solution for A given quantum number Ij 
(and the rest of quantum numbers, which stay the same). This way we can sample Ij from the Pa distribution given 
all other quantum numbers li^j are fixed. Then we iterate the procedure with index j running from 1 to A'' several 
times. Using this algorithm we obtain a set of statistically independent states pooled with probability cx A\ [90| . 

In theory one can do Monte-Carlo summation not only over A states, but over /i and v as well. Nevertheless, 
we found such approach impractical because of a much slower convergence stemmed from very poor cancellation of 
various harmonics. 



C. Results 



1. Ground state 



First we would like to focus on correlation functions for the ground state. We have several reasons for that. First, 
exact analytical solution for such correlation functions is known (j78p , hence we can at least partially verify the validity 
of our approach. Second, we do not have to sum over A and v, so computations are very fast and precise. Therefore 
we can use this case to investigate the effects of our approximations — what finite values of interaction constant c 
can be considered as infinite, how one particle approximation affects results, and what are the finite size effects of our 
calculations. 

First we looked at spatial correlator, calculated at L = 51, N — 51. For the ground state it is translationally 
invariant both in space and time. On Figs. [T]we show a theoretical correlator along with ones calculated at c = 100, 
c — 1000, and c — cx) for one particle processes only. We see that for c = 1000 the correlator is almost indistinguishable 
from the theoretical c — >■ 00 limit, though c = 100 slightly deviates from that limit. Hence all our further c — ^ 00 
results apply to c > 1000 case as well. Also note that theoretical correlation function perfectly overlaps with one 
calculated for c = 00, hence our choice of system size can be considered as a thermodynamic limit. 

As a separate run we calculated the same correlators for the same finite values of c while taking into account all 
two-particle contributions and major three-particle contributions in addition to the original one-particle ones. We 
found that these contributions are negligible. We conclude that it is safe to ignore them in the large c limit. 

On Fig. [5] we show a temporal correlation function for a ground state for c — > 00 as well as a theoretical value given 
by (|78p . We see that they perfectly overlap except for the tiny region around x=0 because of momentum cutoff in 
our numeric calculations. 



2. Delta function in momentum space 



We consider the case when before interactions are switched on all particles are initially in the same state with a 
given momentum p. Such state can be prepared experimentally from the condensate at rest using a Bragg pulse. 

Before we proceed with calculations, we would like to take a closer look at overlap factors A\ for this state (I64p . 
We are starting with the equation for the projection of our state onto eigenfunctions (TG) of the Hamiltonian 
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FIG. 1: Upper panel: calculated and theoretical spatial correlation function for a ground state for c = 100, c — 1000, and 
c — >■ oo. Lower panel: zoom of the upper panel. 



(q)|vl/g,Q^. Ignoring coefficients independent of the TG state and non-singular a,t q — 0, we have for odd |l03l | N 

Note that our system is discretely translationally invariant in k space: we can shift all and q by the same constant 
for arbitrary n, change quantum numbers Ij by n, and all the equations (|45p - (j64p will hold. 

Neglecting finite size effects (resulting in discretization), we can choose n such that q — ^n. Then we shift all Xi 
by q, and consider the overlap ([M)) factor in the g — ?> limit. 

Sinus in the expression (|89|) results in a order zero in the numerator in (7 — >■ oo limit. Zeroes of the denominator 
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FIG. 2: Theoretical and calculated temporal correlator for a ground state in the c — oo limit. Real, imaginary parts, and 
absolute values are shown. Only three curves are visible because theoretical and calculated plots perfectly overlap. 
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FIG. 3: Three-dimensional image (coordinate, time and intensity) of the correlation function for the initial state corresponding 
to the delta peak in momentum space. 



at g = are given by one of A.^ = in the first product, and all the pairs < j such that — —Xj in the 

second product. Maximum order of zero we can achieve in the denominator is when we have one zero Ai, and for each 
of the rest of lambdas we have one exactly opposite to it. For such configuration zero's order of the denominator is 
also ^^^Y^- In this case they cancel resulting in a non-zero contribution. Any other configurations will result in the 
lower order of zero of the denominator, hence zero overlap. We conclude that all the states with non-zero overlap are 
symmetric. One can easily repeat this calculation for even N and show that states in this case are symmetric as well. 
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FIG. 4: Left: spatial correlator for the initial state corresponding to the delta function in momentum space for different to. 
Right: zoom of the left panel. 



We can draw an immediate conclusion about the structure of the sum (ISS]) we use to calculate correlation function. 
Remember that for an infinite c (A|i^|/i) has non-zero elements only if states A and differ by at most one quantum 
number ((73|) . Therefore in the summation ((88| A is symmetric, /i has one quantum number changed relative to A, 
and V, being symmetric, has one quantum number changed relative to ^. But it is possible only if either |A) = \v), 
or if the changed quantum number in |A) — >■ is symmetric to the corresponding change for — >■ \v). I.e. if first 
we change li — >■ than later we should change either 7^' — > /j or — = Ij — > /j = — As a result for a given |A) a 
selection of allows only two possible states for \v), strongly reducing the phase space we need to consider. 

Using the symmetry we can also simplify the expression for overlap factors (j64p . Because in the sum (|88p normaliza- 
tion of does not matter, we drop all the constant terms and can consider ([55)1 instead. Let's denote n — {N — l)/2 
(N is odd), and indexes run from — n to n. Then we can rewrite the numerator as 

n (A,-A,)2(A,+A,f H ~2\] (90) 

0<j<k<n 0<j<n 



and the denominator as 



so that (|M)) simplifies to 



(-9)^2^ n (A.-A.)^A,+A.)^ n (91) 

0<j<fc<ri 0<j<ri 



Ax (X -. (92) 

llAj>0 

In case of interest one can easily restore the constant coefficient and also prove that this expression is valid for even 
N as well. 

Using the conclusions drawn above about relative structure of A, /i, and states, and using expression (|92)) for 
overlap factors, we perform a Monte Carlo summation over states for N — 101, L = 101. On Fig. |4]we show a spatial 
correlation (p(xo, ^0)^(2^0 + x, to)) as a function of x for xq = and different values of Iq. Basically the plot represents 
time evolution of x correlator starting from the moment we momentarily switch on a very strong interaction. We 
can see that at the moment to = the entire correlation function is not distorted except for x = 0. This is because 
when we immediately switch on the interaction, particles stay where they are. They just don't have any time to move 
anywhere, hence the original non-interacting correlation function, which is flat, is preserved (small distortions around 
X = are artifacts of our calculations; if to take into account their errors, they are indistinguishable from 1). But as 
the time goes by, repelling interaction generates a wave, which transfers the matter away from the particle sitting at 
X = 0. On the figure we observe the evolution of this wave. When the time becomes sufficiently large {t > 1), the 
correlation function stabilizes and does not evolve anymore. 

We also calculate a temporal correlator {p{xo,to)p{xo,to + t)). Figure [5] shows real and imaginary parts of the 
correlator. Because the state is translationally invariant, correlator does not depend on xq. Interestingly, it also 
almost does not depend on Iq. All the plots for different to overlapped, and the difference was invisible. On this plot 
we also observe two "effects" , which are artifacts of our calculation scheme because of introduced momentum cutoffs: 
oscillations in the imaginary part and real part not going exactly to zero at i = 0. 
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FIG. 5: Real and imaginary parts of a temporal correlator for the initial state corresponding to the delta function in momentum 
space. 



Now we proceed with the initial state — a Gaussian pulse in space, and hence in momentum space. Overlap factors 
for such state are given by the expression ([68]) . We neglect finite size effects by considering the hmit L ^ cr, for which 
expression up to a constant factor, can be written as 



In this case we do not have a nice symmetry we had in the previous section, so we have to sum over all one particle 
II states and two particle v states using expression ([55)1 . Summation over all two particle processes substantially 
increases the computational time, so in order to partially compensate for it we used smaller values for problem 
parameters L = 21 and N = 21. For a packet width parameter we use value a — 15/(27r). 

On Fig.[7]we present time evolution of spatial correlator {p{xQ,tQ)p{xQ +x,t())) for xq = and different fixed times 
to. We also show errors of our calculations to distinguish true correlation function features from noise. 

We would like to emphasize several features visible on these plots. First, for to = 0, as expected, we see a dip 
at a; = 0. Interestingly, the width of this dip is finite and is not determined by the momentum cutoff. From 
computational point of view we can explain it by the fact that the overlap factor (j93p for Gaussian state strongly 
inhibits high harmonics, therefore we do not have momenta high enough to make this dip infinitesimally narrow. This 
is a finite size effect because of a small number of particles. As parameters L, N, and a increase proportionally, width 
of the dip goes to zero. Second, oscillations of the correlator at = around the dip a:o = are not determined by 
the cutoff and much bigger than the error of our calculations. This is again the effect of a finite size system. Third, 
as we've already seen in previous section, for times to ^ 1 the correlation function stabilizes. But the process of 
stabilization is quite remarkable. We do not have a widening of a Gaussian peak. Instead it looks like the peak has 
the same width, but its height decreases with time. Probably, qualitative explanation lies in the fact that high density 
regions have high interaction energies, hence in these regions we have fast particles which quickly leave the region. 
I.e. we have a ballistic scenario, which differs from the diffusion scenario for which we would expect a widening of a 
Gaussian peak with time. 

Unfortunately, due to some peculiarities of computations we had to conduct, we were unable to produce reliable 
results for a temporal correlator. Results, though reproducing an overall shape of the correlator, were very noisy 
because of beats attributed to non-exact cancellation of various harmonics. We do not provide a plot for these results. 



3. Gaussian pulse 




(93) 




FIG. 7: Spatial correlator for different to for a Gaussian initial state and its error. Two lines for each value of to lie one standard 
deviation above and below relative to the average value. 
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V. DISCUSSION AND SUMMARY 



In this paper we addressed a problem of non-equilibrium time evolution of one dimensional Bose gas with contact 
interaction from the general perspective of dynamics of integrable systems. Several approaches have been proposed 
recently for analyzing time evolution of integrable many-body models, including Quantum Inverse Scattering method, 
the formalism of Intertwining Operator, and Extended Conformal Symmetry. After critically reviewing these ap- 
proaches we concentrated on conceptually the simplest method, based on using Bethe ansatz to decompose initial 
states into precise eigenstates of the interacting Hamiltonian using the intertwining operator and then using form 
factor expansion to calculate time evolution of correlation functions. The main difficulty of this method is that it 
requires summation over a large number of intermediate states. In this paper we focused on the regime of strong 
repulsive interactions between bosons and developed an efficient numerical procedure for performing summation over 
intermediate states. We analyzed two types of initial states: all particles having zero momentum and all particles in 
a gaussian wavepacket in real space. In both cases we find a non-trivial time evolution of the second order coherence 
g2{xi, X2,t) = {p{xi,t)^ p{x2,t)), which reflects the intrinsic dichotomy between initial states and the Hamiltonian. 
Initial condensate states exhibit bunching at short distances, whereas strong repulsion in the Hamiltonian introduces 
strong antibunching. 

For the initial state which has all particles in a state with zero momentum, results for (72 are summarized in Fig. [31 
Antibunching at shortest distances is present at all times reflecting repulsion between particles. A more striking 
feature is the appearance of bunching and oscillations in g2ixi — X2,t) at intermediate time and length-scales. At 
transient times 52 has Friedel like oscillations as a function of X1 — X2, which disappear at longer times. Crystallization 
of the TG gas in the process of non-equilibrium time evolution was also discussed in Ref 20]. At longer times the form 
of 52 is qualitatively similar to what one flnds for the equilibrium Lieb Liniger model with intermediate interactions 
at high temperatures Fig. [3] also provides some support to the idea of light cone formation discussed in Ref. [9l| . 
Some of these features are similar to quench dynamics in other integrable systems [s^, • 

For the initial state which has all particles in a gaussian wavepacket in real space, time evolution of 52 is shown in 
Figs. I6I7I As in the previous case we observe antibunching at the shortest distances. For intermediate time scales it 
is followed by a pronounced peak in 172 reflecting dynamic bunching. This system would be naturally characterized 
by the time dependent short range correlation length, which increases with time, reflecting expansion of the system in 
real space. Oscillations in 52 which we observe in this case are small and are most likely related to the finite number 
of particles used in the analysis. Hence transient states of this system would be more natural characterized as a liquid 
rather than a crystal. 

Predictions made in our paper for the behavior of 52 should be possible to test in experiments with ultracold atoms 
and photons in strongly non-linear medium. 
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VII. APPENDIX A: INVERSE SCATTERING TRANSFORM AND ALGEBRAIC BETHE ANSATZ - 

BASIC CONCEPTS 



The nonlinear Schrodinger Hamiltonian represents the simplest example of a system solvable by the algebraic Bethe 
ansatz [ss', . For review and many details see 4] . Here we overview the basic construction with connection to the 
inverse scattering transform in order to fix notations and for the sake of self-consistency. 

The Zakharov-Shabat method starts from transforming the field 'i>{x,t = 0) of nonlinear Hamiltonian (HI) at time 
t ^ into a set of " scattering data" given by the following linear problem 

i^^x,o = Qix.OHx.O (94) 



where 
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The solution ^) of this equation is defined by the condition | ^' (x) | — )■ as x — ±cxd and by the properties of Jost 
solutions, 

(97) 



^(^)g.5./2 j ^.^-oo i^;^2(x,e);™ Vie"'*"/' 

Rewriting the linear equation above with the boundary conditions as an integral equation gives the Gelfand-Levitan- 
Marchenko equation 

/oo 
dyd{y - x)e^^y^\y)xi{y. i)e''^^l\ (98) 
-OO 

/oo 
dy0(2/-a;)e-'«^Xi(2/,C)*(y)e'^«/' (99) 
-OO 

which can be solved iteratively in powers of ^fc 

oo n ^ n 

Xi(x,Oe"*/' = c" n / n ^2^J^(^o >yi> ■■■yn>xn> x)e*«(^"=o^'-^"=i (100) 

n=0 'i=0 j = l 

X *t(2.^)...^t(3,^)^(yj...^(y^)^ (101) 

oo n— 1 ^ n 

X2(x, Oe"'"*/' = 1 + E ^" n / ^^2;, n '^yj^C^^o > yi > • • • x„_i > y„ > a;)e^«(^-«' "--^"-^ (102) 

X ^^xo)---^Hxr,^i)^i'{y„)---^i'{yi) (103) 

where 0(xi > X2 > • • ■ > x„) = 0(a;i — X2)0{x2 — x^) ■ ■ ■ 9{xn-i — Xn). This solution can be written as a solution for 
the scattering data, A{^) and B{^). Finally, defining the reflection operator as 

RiO^mOr'BiO (104) 

one obtains the expansion Operators R{^) satisfy the Zamolodchikov-Faddeev algebra In this formalism, the 
Bethe ansatz states are constructed by application of Zamolodchikov-Faddeev operators to the pseudovacuum state 

|0), 

\^N)BA=RH^l).--RH^Nm (105) 

These states are complete for the case of repulsion. Using the expansion ([2]) one can show [32! | that the special ordered 
initial state of the form 

\Mord = 0{xi >X2> ...> XN)\'i\xi)^\x2) . . . *^(xjv)|0) (106) 

is equal to the following state 

R\xi)R\x2)...R\xn)\0) (107) 

where 

Rix) = / §e"«i?(e). (108) 

We note that the evolution of this state can be obtained explicitly, as discussed in Section II. The case of attraction 
is discussed in Appendix C. 

The transfer matrix is a central ingredient of the construction of the inverse scattering transform on the finite 
interval on the lattice (one performs a space discretization procedure for the LL model first). It can be represented 
as a matrix 



(109) 
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The transfer matrix is constructed as a product of monodromy matrices ij(^) for each site j, T(^) = Y[j which 
in the case of NS problem has the form of matrix Q which is one of matrix forming the Lax pair. Trace of the transfer 
matrix commute for different ^ and is a generator of integrals of motion. The integrability condition is expressed by 
the special property 

- am) ® TiO) = ina ® t{o)r{^ - eo (no) 

where the matrix R{^, is a solution of the Yang-Baxter equation, which for the case of NS problem belongs to the 
class of 6-vertex model 

RuiORisi^ + m)^23(m) = R23{^^)Rl3{^ + ^^)Rl2{0 (m) 

Rab=l3Iab + aPab, a ^ ^\ ^\ , (3 = . (112) 

Ca - ?6 - «C ?a - ?6 - «C 

where lat and Pab arc identity and permutation operators acting on the tensor product of two single-particle spaces 
indexed by a and b. In the finite size lattice formalism the diagonal entries of the matrix (more precisely its trace 
'''{0 = TrT{^) — A{^) + D{^)) generate a conserved quantities of the model whereas the off-diagonal elements B{^) 
and C(^) act as creation and annihilation operators of pseudoparticles. The n-particle Bethe ansatz eigenvectors are 
constructed using these off-diagonal elements of the transfer matrix as 

|*„({A„})) - B{X,)B{X2) . . . B(A„)|0) (113) 

whereas the bra- vectors are constructed as 

(*„({A„})| = (0|C(Ai)C(A2) . . . C(A„) (114) 

The algebraic Bethe ansatz deals with diagonalization of the trace r(f) of the transfer matrix, thus solving the 
eigenvalue equation t{X)\BA) = A(A, {Xa})\BA) . This can be done using the commutation relations between operators 
A{X), B{X),C{X), D{X) which can be deduced from the Yang-Baxter algebra. This procedure leads to expression for 
the A(A, {Xa}) = a(A) 0"=! /(A ^ Xa) + /3(A) Ha^i /(Aa — A) which can be further related to the energy and to the 
consistency relation, known as a Bethe ansatz equations for the momenta Xa of the pseudo particles, 

"(Aa) _ TT Afc - Aa 



where the function /(A) is determined by the elements of the i?-matrix. The connection with the continuum version 
of the NLS model is established by a limiting procedure of lattice size going to zero while keeping the density constant 
on a finite interval. _ 

Due to the work of Sklyanin 94] it is possible to formulate even more general NS problem which includes the 
boundaries. The boundary problems can be divided into two categories, soliton-preserving and soliton-non-preserving. 
For recent progress in solution of the boundary NS model see [95| . 
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